2023 年 数维杯（A题）国际大学生数学建模挑战赛 建模秘籍&文章代码思路大全

铛铛！小秘籍来咯！
小秘籍希望大家都能轻松建模呀，数维杯也会持续给大家放送思路滴~
抓紧小秘籍，我们出发吧~
来看看 数维杯国际大学生数学建模挑战赛的A题！
 

铛铛！小秘籍来咯！ 小秘籍希望大家都能轻松建模呀，数维杯也会持续给大家放送思路滴~ 抓紧小秘籍，我们出发吧~ 来看看 数维杯国际大学生数学建模挑战赛的A题！
 
## 问题重述
 
1、利用给定参数，建立复合直升机的俯仰矩表达式，并发展俯仰角变化的模型。给出在初始飞行条件下（初始飞行高度3000米，前进速度80 m/s，垂直上升速度2 m/s，初始俯仰角为0度等参数）飞行机体在5秒、10秒和20秒时的姿态角度。
 
2、基于题目给定的参数，建立复合直升机的横滚、俯仰和偏航力矩的表达式，并发展姿态角度变化的模型。给出在初始飞行条件下（初始飞行高度3000米，前进速度80 m/s，垂直上升速度0.2 m/s，初始横滚角、俯仰角和偏航角等参数）飞行机体在5秒、10秒和20秒时的姿态角度。
 
3、基于直升机在低速和高速飞行中的机动特性，设计每个组件的机动幅度，以满足直升机在水平飞行任务（零姿态角）下的需求。考虑两种情况：初始飞行高度3000米，前进速度80 m/s，垂直上升速度0.2 m/s，以及初始飞行高度3000米，前进速度180 m/s，垂直上升速度0.2 m/s。
 
4、在直升机进行加速机动任务的情况下，飞行速度的前向分量在20秒内均匀增加从80 m/s到180 m/s（假设加速度完全由螺旋桨推进器提供），请设计每个控制输入的动态值，以实现前向加速度和水平飞行（零姿态角）。考虑直升机在低速和高速飞行中的机动特性，初始飞行高度为3000米，垂直上升速度为0.2 m/s。
### 问题一：详细建模思路
 
#### 1. **Pitch Moment 表达式的建立：**
 
考虑俯仰矩 $M_{pitch}$，它是直升机俯仰运动时产生的力矩。这个力矩来自多个组件：
 
- **Coaxial Rigid Rotors:** 这部分力矩可以通过以下方式建模：
 $M_{rotor} = K_r \cdot \rho \cdot A \cdot V^2 \cdot R \cdot \Omega^2$
 其中，$K_r$ 是 rotor moment factor，$\rho$ 是空气密度，$A$ 是转子面积，$V$ 是飞行速度，$R$ 是转子旋翼半径，$\Omega$ 是转子角速度。
 
- **Propeller Thrusters:** 螺旋桨推进器引起的力矩可表示为：
 $M_{prop} = \text{（根据螺旋桨工作状态和转速计算的力矩）}$
 
- **Horizontal Tails:** 水平尾翼通过改变相对角度引起的力矩可表示为：
 $M_{h} = C_{m_h} \cdot q \cdot S_h \cdot d_{h}$
 其中，$C_{m_h}$ 是水平尾翼力矩系数，$q$ 是动压，$S_h$ 是水平尾翼面积，$d_{h}$ 是水平尾翼到重心的距离。
 
- **Vertical Tails:** 垂直尾翼引起的力矩可表示为：
 $M_{v} = C_{m_v} \cdot q \cdot S_v \cdot d_{v}$
 其中，$C_{m_v}$ 是垂直尾翼力矩系数，$S_v$ 是垂直尾翼面积，$d_{v}$ 是垂直尾翼到重心的距离。
 
综合考虑以上组件的力矩，可以得到总的俯仰矩表达式：
$M_{pitch} = M_{rotor} + M_{prop} + M_{h} + M_{v}$
 
#### 2. **牛顿-欧拉方程的建立：**
 
利用牛顿-欧拉方程，可以将直升机的运动和外部力矩联系起来。一般形式如下：
 
$I \cdot \dot{\omega} = \tau - \omega \times (I \cdot \omega)$
 
其中，$I$ 是惯性矩阵，$\dot{\omega}$ 是角速度的导数，$\tau$ 是总外部力矩，$\omega$ 是角速度。这个方程需要根据直升机的几何结构和组件的动力学影响进行详细的分析和建模。
 
#### 3. **参数表示：**
 
将各组件的影响表示成符号变量，例如：
 
- $K_r$：共轴刚性转子的力矩系数
- $C_{m_h}$：水平尾翼的力矩系数
- $C_{m_v}$：垂直尾翼的力矩系数
 
这样可以更清晰地表示模型中的各项参数。
 
#### 4. **飞行模式切换逻辑：**
 
建立逻辑以在低速和高速两种模式之间切换。可以考虑使用阈值条件，如飞行速度 $V$ 是否超过某个值。
 
#### 5. **数学模型整合：**
 
将建立的各个方程整合为一个完整的数学模型。考虑牛顿-欧拉方程、俯仰矩方程和飞行模式切换逻辑。
 
#### 6. **数值模拟验证：**
 
使用数值模拟工具，例如Matlab或Python中的NumPy，对建立的数学模型进行验证。通过模拟不同飞行条件下的俯仰运动，检查模型的合理性和稳定性。
 
```python
import numpy as np
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt
 
# 定义直升机动力学方程
def helicopter_dynamics(t, y, rho, A, R, Omega, C_m_h, C_m_v, S_h, S_v, d_h, d_v, K_r):
   phi, theta, psi, p, q, r = y
 
   # 计算俯仰矩
   M_pitch = K_r * rho * A * R**5 * Omega**2
 
   # 计算水平尾翼和垂直尾翼的力矩
   M_h = C_m_h * 0.5 * rho * q**2 * S_h * d_h
   M_v = C_m_v * 0.5 * rho * r**2 * S_v * d_v
 
   # 牛顿-欧拉方程
   I_x = 100.0 # X轴惯性矩
   I_y = 150.0 # Y轴惯性矩
   I_z = 200.0 # Z轴惯性矩
 
   p_dot = (I_y - I_z) * q * r / I_x - M_pitch / I_x
   q_dot = (I_z - I_x) * p * r / I_y + M_h / I_y
   r_dot = (I_x - I_y) * p * q / I_z + M_v / I_z
 
   # 返回角速度的导数
   return [p_dot, q_dot, r_dot, 0, 0, 0]
 
# 定义初始条件
initial_conditions = [0.0, 0.0, 0.0, 0.1, 0.1, 0.1]
 
# 定义时间范围
time_span = (0, 10)
 
# 定义其他参数
rho = 1.225 # 空气密度
A = 10.0 # 转子面积
R = 5.0 # 转子半径
Omega = 200.0 # 转子角速度
C_m_h = 0.1 # 水平尾翼力矩系数
C_m_v = 0.05 # 垂直尾翼力矩系数
S_h = 2.0 # 水平尾翼面积
S_v = 1.0 # 垂直尾翼面积
d_h = 3.0 # 水平尾翼到重心的距离
d_v = 2.0 # 垂直尾翼到重心的距离
K_r = 0.02 # 共轴刚性转子的力矩系数
 
# 数值模拟
solution = solve_ivp(
   fun=lambda t, y: helicopter_dynamics(t, y, rho, A, R, Omega, C_m_h, C_m_v, S_h, S_v, d_h, d_v, K_r),
   t_span=time_span,
   y0=initial_conditions,
   method='RK45',
   dense_output=True
)
 
# 获取结果
t_eval = np.linspace(time_span[0], time_span[1], 1000)
result = solution.sol(t_eval)
 
# 绘制结果
plt.figure(figsize=(10, 6))
```
### 问题二：
#### 分析问题：
 
1. **Roll, Pitch, Yaw Moments 表达式：**
  - 分别建立直升机在给定飞行条件下的滚转矩、俯仰矩和偏航矩的表达式。
  - 使用力矩平衡方程，考虑共轴刚性转子、螺旋桨推进器、升降舵和方向舵对各个力矩的贡献。
 
2. **动力学方程：**
  - 利用牛顿-欧拉方程建立直升机的动力学方程，同时考虑在不同飞行模式下的动力学行为。
 
3. **参数表示：**
  - 用符号表示每个组件的影响，如 rotor moment factor、horizontal tail moment coefficient、vertical tail moment coefficient 等。
 
4. **不同飞行模式的切换：**
  - 考虑低速和高速两种飞行模式，根据飞行速度的变化切换相应的动力学模型。
 
5. **力矩和姿态角的关系：**
  - 将建立的力矩方程与姿态角（滚转、俯仰、偏航）的变化关系联系起来。
 
#### 具体建模步骤：
 
1. **Roll Moment 表达式的建立：**
  - 定义滚转矩 $M_{roll}$，考虑共轴刚性转子、螺旋桨推进器、升降舵和方向舵的力矩。
  - 使用力矩平衡方程表达 $M_{roll}$。
 
2. **Pitch Moment 表达式的建立：**
  - 与问题一类似，建立俯仰矩 $M_{pitch}$ 的表达式。
 
3. **Yaw Moment 表达式的建立：**
  - 定义偏航矩 $M_{yaw}$，考虑共轴刚性转子、螺旋桨推进器、升降舵和方向舵的力矩。
  - 使用力矩平衡方程表达 $M_{yaw}$。
 
4. **牛顿-欧拉方程的建立：**
  - 使用牛顿-欧拉方程，将飞行器的姿态角和角速度与外部力和力矩联系起来。
  - 根据直升机的几何结构和组件的动力学影响，建立不同飞行模式下的牛顿-欧拉方程。
 
5. **参数化表示：**
  - 将涉及的参数进行符号表示，例如 rotor moment factor ($K_r$)、horizontal tail moment coefficient ($C_{m_h}$)、vertical tail moment coefficient ($C_{m_v}$) 等。
 
6. **飞行模式切换逻辑：**
  - 根据飞行速度的阈值，建立逻辑以在低速和高速两种模式之间切换。
 
7. **数学模型整合：**
  - 将上述步骤中建立的方程整合成一个综合的数学模型，描述直升机在不同飞行条件下的滚转、俯仰和偏航力矩以及姿态角的变化。
 
8. **数值模拟验证：**
  - 使用数值模拟工具，验证建立的数学模型在给定飞行条件下的准确性和稳定性。
 
```python
import numpy as np
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt
 
# 定义直升机动力学方程
def helicopter_dynamics(t, y, rho, A, R, Omega, C_m_h, C_m_v, S_h, S_v, d_h, d_v, K_r):
   phi, theta, psi, p, q, r = y
 
   # 计算滚转、俯仰、偏航力矩
   M_roll = K_r * rho * A * R**5 * Omega**2 # 仅作示例，实际需要更详细的建模
   M_pitch = C_m_h * 0.5 * rho * q**2 * S_h * d_h # 仅作示例，实际需要更详细的建模
   M_yaw = C_m_v * 0.5 * rho * r**2 * S_v * d_v # 仅作示例，实际需要更详细的建模
 
   # 牛顿-欧拉方程
   I_x = 100.0 # X轴惯性矩
   I_y = 150.0 # Y轴惯性矩
   I_z = 200.0 # Z轴惯性矩
 
   p_dot = (I_y - I_z) * q * r / I_x - M_roll / I_x
   q_dot = (I_z - I_x) * p * r / I_y - M_pitch / I_y
   r_dot = (I_x - I_y) * p * q / I_z - M_yaw / I_z
 
   # 返回角速度的导数
   return [p_dot, q_dot, r_dot, 0, 0, 0]
 
# 定义初始条件
initial_conditions = [0.0, 0.0, 0.0, 0.1, 0.1, 0.1]
 
# 定义时间范围
time_span = (0, 10)
 
# 定义其他参数
rho = 1.225 # 空气密度
A = 10.0 # 转子面积
R = 5.0 # 转子半径
Omega = 200.0 # 转子角速度
C_m_h = 0.1 # 水平尾翼力矩系数
C_m_v = 0.05 # 垂直尾翼力矩系数
S_h = 2.0 # 水平尾翼面积
S_v = 1.0 # 垂直尾翼面积
d_h = 3.0 # 水平尾翼到重心的距离
d_v = 2.0 # 垂直尾翼到重心的距离
K_r = 0.02 # 共轴刚性转子的力矩系数
 
# 数值模拟
solution = solve_ivp(
   fun=lambda t, y: helicopter_dynamics(t, y, rho, A, R, Omega, C_m_h, C_m_v, S_h, S_v, d_h, d_v, K_r),
   t_span=time_span,
   y0=initial_conditions,
   method='RK45',
   dense_output=True
)
 
# 获取结果
t_eval = np.linspace(time_span[0], time_span[1], 1000)
result = solution.sol(t_eval)
 
# 绘制结果
plt.figure(figsize=(12, 8))
plt.subplot(3, 1, 1)
plt.plot(t_eval, result[0], label='Roll Angle (phi)')
plt.legend()
plt.subplot(3, 1, 2)
plt.plot(t_eval, result[1], label='Pitch Angle (theta)')
#见完整代码
```
### 问题三：建模思路
#### 问题分析：
 
1. **低速飞行模式：**
  - 根据题目描述，低速模式下主要由共轴刚性转子和螺旋桨推进器控制。
  - 建立在低速模式下保持零姿态角的条件下，各组件的操纵量。
 
2. **高速飞行模式：**
  - 在高速模式下，主要由螺旋桨推进器、升降舵和方向舵控制。
  - 建立在高速模式下保持零姿态角的条件下，各组件的操纵量。
 
3. **操纵量与姿态角关系：**
  - 将各组件的操纵量与姿态角的关系建立起来，考虑其对滚转、俯仰和偏航的影响。
 
4. **零姿态角条件下的操纵量设计：**
  - 在两种飞行模式下，设计使得直升机保持零姿态角所需的操纵量。
 
5. **数学模型整合：**
  - 将上述步骤中建立的关系整合成一个综合的数学模型，描述直升机在不同飞行模式下实现零姿态角的操纵量。
 
6. **数值模拟验证：**
  - 使用数值模拟工具，验证建立的数学模型在给定飞行条件下的准确性和稳定性。
 
#### 具体建模步骤：
 
1. **低速飞行模式的建模：**
  - 考虑共轴刚性转子和螺旋桨推进器对低速模式下姿态角的影响。
  - 设计合适的操纵量，使得在低速飞行时保持零姿态角。
 
2. **高速飞行模式的建模：**
  - 考虑螺旋桨推进器、升降舵和方向舵对高速模式下姿态角的影响。
  - 设计合适的操纵量，使得在高速飞行时保持零姿态角。
 
3. **操纵量与姿态角的关系建立：**
  - 将操纵量与姿态角的关系表达成数学方程。
  - 考虑不同操纵量对姿态角的影响，如俯仰操纵量对俯仰角的影响等。
 
4. **零姿态角条件下的操纵量设计：**
  - 对于低速和高速模式，设计操纵量的变化规律，以保持零姿态角。
  - 可能需要进行迭代和调整，确保在不同飞行模式下都能有效地保持零姿态角。
 
5. **数学模型整合：**
  - 将两种飞行模式下的建模结果整合，形成直升机在不同飞行条件下的操纵量与姿态角关系的综合模型。
 
6. **数值模拟验证：**
  - 使用数值模拟工具，验证建立的数学模型在给定飞行条件下的准确性和稳定性。
  - 考虑不同初始条件和飞行任务，以确保模型的通用性。
 
```python
import numpy as np
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt
 
# 定义直升机动力学方程，包括低速和高速两种模式
def helicopter_dynamics(t, y, mode, control_inputs):
   phi, theta, psi, p, q, r = y
 
   # 根据飞行模式选择不同的操纵量
   if mode == 'low_speed':
       elevator, rudder = control_inputs['low_speed']
   elif mode == 'high_speed':
       elevator, rudder = control_inputs['high_speed']
   else:
       raise ValueError("Invalid flight mode")
 
   # 按照题目中的公式编写滚转、俯仰和偏航的动力学方程
   # ...
 
   # 返回角速度的导数
   return [p_dot, q_dot, r_dot, 0, 0, 0]
 
# 定义初始条件
initial_conditions = [0.0, 0.0, 0.0, 0.1, 0.1, 0.1]
 
# 定义时间范围
time_span = (0, 10)
 
# 定义操纵量
control_inputs = {
   'low_speed': [-0.01, 0.02], # 低速模式下的操纵量
   'high_speed': [0.01, -0.02] # 高速模式下的操纵量
}
 
# 数值模拟，分别模拟低速和高速两种飞行模式
solution_low_speed = solve_ivp(
   fun=lambda t, y: helicopter_dynamics(t, y, 'low_speed', control_inputs),
   t_span=time_span,
   y0=initial_conditions,
   method='RK45',
   dense_output=True
)
 
solution_high_speed = solve_ivp(
   fun=lambda t, y: helicopter_dynamics(t, y, 'high_speed', control_inputs),
   t_span=time_span,
   y0=initial_conditions,
   method='RK45',
   dense_output=True
)
 
# 获取结果
t_eval = np.linspace(time_span[0], time_span[1], 1000)
result_low_speed = solution_low_speed.sol(t_eval)
result_high_speed = solution_high_speed.sol(t_eval)
 
# 绘制结果
plt.figure(figsize=(12, 8))
plt.subplot(3, 2, 1)
plt.plot(t_eval, result_low_speed[0], label='Roll Angle (phi) - Low Speed')
plt.legend()
plt.subplot(3, 2, 2)
plt.plot(t_eval, result_high_speed[0], label='Roll Angle (phi) - High Speed')
plt.legend()
plt.subplot(3, 2, 3)
plt.plot(t_eval, result_low_speed[1], label='Pitch Angle (theta) - Low Speed')
plt.legend()
plt.subplot(3, 2, 4)
plt.plot(t_eval, result_high_speed[1], label='Pitch Angle (theta) - High Speed')
```
 
 
## **问题四：**
 
在20秒内，直升机的前向飞行速度从80 m/s均匀增加到180 m/s，初始飞行高度为3000米，垂直上升速度为0.2 m/s。要求设计每个控制输入的动态值，以实现前向加速度和水平飞行（零姿态角）。考虑直升机在低速和高速飞行中的机动特性。
 
**建模思路：**
 
1. **机动特性建模：**
  - 在低速和高速模式下，设计共轴刚性转子、螺旋桨推进器、升降舵、方向舵的机动特性。
  - 将机动特性与操纵量之间的关系用数学模型表示。
 
2. **前向加速度的要求：**
  - 根据速度增加的要求，设计一个函数描述速度随时间的变化，计算加速度。
 
3. **控制输入的动态值设计：**
  - 使用机动特性的数学模型，结合计算得到的加速度，设计螺旋桨推进器、升降舵、方向舵的动态值。
  - 注意在整个时间范围内均匀增加速度的要求。
 
4. **数学模型整合：**
  - 将机动特性的数学模型与控制输入的动态值设计整合，形成一个综合的数学模型。
  - 该模型描述了在加速机动任务中，直升机通过调整控制输入，实现前向加速度和水平飞行（零姿态角）的过程。
 
5. **数值模拟验证：**
  - 使用数值模拟工具，验证建立的数学模型在给定飞行条件下的准确性和稳定性。
  - 考虑不同初始条件、加速度任务和飞行要求，以确保模型的通用性。
 
 
```python
import numpy as np
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt
 
# 定义直升机动力学方程，包括低速和高速模式
def helicopter_dynamics(t, y, mode, acceleration_profile):
   # y: [phi, theta, psi, p, q, r]
   phi, theta, psi, p, q, r = y
   
   # 根据问题四的要求，设计螺旋桨推进器的动态值
   thrust = design_thrust(t, acceleration_profile)
   elevator, rudder = design_controls(t, acceleration_profile)
   
   # TODO: 编写直升机动力学方程，考虑低速和高速模式下的影响
   # ...
 
   return [p_dot, q_dot, r_dot, 0, 0, 0]
 
# 设计螺旋桨推进器的动态值，以实现前向加速度
def design_thrust(t, acceleration_profile):
   # TODO: 根据时间和加速度曲线设计螺旋桨推进器的动态值
   # ...
 
# 设计升降舵和方向舵的动态值
def design_controls(t, acceleration_profile):
   # TODO: 根据时间和加速度曲线设计升降舵和方向舵的动态值
   # ...
 
# 定义初始条件
initial_conditions = [0.0, 0.0, 0.0, 0.1, 0.1, 0.1]
 
# 定义时间范围
time_span = (0, 20)
 
# 定义加速度的变化规律
def acceleration_profile(t):
   # TODO: 根据时间设计加速度的变化规律
   # ...
 
# 数值模拟
solution = solve_ivp(
   fun=lambda t, y: helicopter_dynamics(t, y, 'low_speed', acceleration_profile),
   t_span=time_span,
   y0=initial_conditions,
   method='RK45',
   dense_output=True
)
 
# 获取结果
t_eval = np.linspace(time_span[0], time_span[1], 1000)
result = solution.sol(t_eval)
 
# 绘制结果
plt.figure(figsize=(12, 8))
plt.subplot(3, 2, 1)
plt.plot(t_eval, result[0], label='Roll Angle (phi)')
plt.legend()
plt.subplot(3, 2, 3)
plt.plot(t_eval, result[1], label='Pitch Angle (theta)')
```