2023 第九届数维杯A题 建模解析，小鹿学长带队指引全代码文章与思路

我是小鹿学长，就读于上海交通大学，截至目前已经帮200+人完成了建模与思路的构建的处理了～

这回带大家体验一下第九届数维杯A题呀！

问题重述

问题一：建立姿态角模型

背景：

直升机的飞行受到多种力的影响，包括共轴刚性转子、螺旋桨推进器、升降舵、方向舵等。在不同飞行模式下，这些组件对直升机的姿态角产生不同影响。

任务：

1. 建模任务： 基于给定参数，建立直升机的俯仰角模型，考虑共轴刚性转子、螺旋桨推进器、升降舵和方向舵对俯仰角的影响。
2. 模型验证： 提供直升机在初始飞行条件下（海拔3000米，前进速度80 m/s，垂直上升速度2 m/s，初始俯仰、横摆、偏航角分别为0度）的姿态角在5秒、10秒和20秒时的数值。

---

问题二：建立姿态角变化模型

背景：

直升机在飞行过程中，其姿态角会随时间变化。不同组件的力矩会导致滚转、俯仰和偏航的变化。

任务：

1. 建模任务： 基于给定参数，建立直升机滚转、俯仰和偏航角的变化模型，考虑共轴刚性转子、螺旋桨推进器、升降舵和方向舵的影响。
2. 模型验证： 提供直升机在初始飞行条件下（海拔3000米，前进速度80 m/s，垂直上升速度0.2 m/s，初始滚转、俯仰、偏航角分别为0度）的姿态角在5秒、10秒和20秒时的数值。

---

问题三：设计机动特性

背景：

直升机在低速和高速飞行模式下，需要满足水平飞行任务，即保持零姿态角。各组件的机动特性在不同飞行模式下起着关键作用。

任务：

1. 机动设计： 根据低速和高速飞行特性，设计共轴刚性转子、螺旋桨推进器、升降舵和方向舵的机动幅度，以满足水平飞行任务。
2. 模型验证： 提供直升机在两种飞行模式下（初始飞行条件相同）的姿态角在5秒、10秒和20秒时的数值。

---

问题四：加速机动任务设计

背景：

直升机需要进行加速机动任务，通过调整控制输入，在20秒内将前向飞行速度从80 m/s均匀增加到180 m/s。

任务：

1. 动态设计： 设计每个控制输入的动态值，以实现前向加速度和水平飞行（零姿态角），考虑低速和高速飞行中的机动特性。
2. 模型验证： 提供直升机在初始飞行条件下（海拔3000米，垂直上升速度0.2 m/s）的姿态角在5秒、10秒和20秒时的数值。

---

建模思路

问题一

1. 基本动力学方程：

在直升机的俯仰运动中，我们可以使用牛顿-欧拉方程来表示：

$I_y\cdot \dot{q}=\sum {\tau }_y$

其中：

• $I_y$ 是绕俯仰轴的惯性矩；
• $\dot{q}$ 是俯仰角速度；
  -$\sum {\tau }_y$ 是所有在俯仰方向上的力矩之和。

2. 各组件的力矩：

考虑到直升机的各个组件对俯仰角的影响，我们可以将力矩项扩展为各组件的贡献之和：

$\sum {\tau }_y={\tau }_{\text{rotors}}+{\tau }_{\text{propeller}}+{\tau }_{\text{elevator}}+{\tau }_{\text{rudder}}$

其中：
-${\tau }_{\text{rotors}}$ 是由共轴刚性转子产生的气动力矩；
-${\tau }_{\text{propeller}}$是由螺旋桨推进器产生的推力力矩；

• ${\tau }_{\text{elevator}}$ 是由升降舵产生的控制力矩；
• ${\tau }_{\text{rudder}}$是由方向舵产生的控制力矩。

3. 具体建模：

具体建模时，需要使用气动和动力学的原理，将各组件的力矩进行数学表达。这可能会包括考虑旋翼气动力矩、螺旋桨推进器的推力、升降舵和方向舵的控制输入等。

例如，对于共轴刚性转子的气动力矩，可以考虑以下一般形式：

${\tau }_{\text{rotors}}=k_{\text{rotors}}\cdot \dot{q}$

其中，$k_{\text{rotors}}$ 是与气动特性相关的比例常数。

类似地，对于其他组件的力矩也需要进行具体建模。

4. 数学模型整合：

将以上的各组件力矩的具体表达式代入基本动力学方程，得到完整的俯仰角动力学方程。

$I_y\cdot \dot{q}={\tau }_{\text{rotors}}+{\tau }_{\text{propeller}}+{\tau }_{\text{elevator}}+{\tau }_{\text{rudder}}$

5. 数值模拟验证：

使用数值模拟工具，如数值积分器，对建立的数学模型进行验证，观察在给定初始条件下，直升机的俯仰角随时间的变化。

import numpy as np
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt

# 定义直升机俯仰角动力学方程
def pitch_dynamics(t, y):
    # y: [pitch, pitch_rate]
    pitch, pitch_rate = y
    
    # 参数设定（示例值，实际应根据直升机具体参数进行调整）
    I_y = 1000.0  # 绕俯仰轴的惯性矩
    k_rotors = 10.0  # 共轴刚性转子的气动力矩比例常数
    k_propeller = 5.0  # 螺旋桨推进器的推力力矩比例常数
    k_elevator = 2.0  # 升降舵的控制力矩比例常数
    k_rudder = 2.0  # 方向舵的控制力矩比例常数
    
    # 各组件的力矩
    rotors_torque = k_rotors * pitch_rate
    propeller_torque = k_propeller * pitch_rate
    elevator_torque = k_elevator * pitch_rate
    rudder_torque = k_rudder * pitch_rate
    
    # 动力学方程
    pitch_acceleration = (rotors_torque + propeller_torque + elevator_torque + rudder_torque) / I_y

    return [pitch_rate, pitch_acceleration]

# 初始条件
initial_conditions = [0.0, 0.0]  # 初始俯仰角为0，初始角速度为0

# 时间范围
time_span = (0, 20)

# 数值模拟
solution = solve_ivp(
    fun=pitch_dynamics,
    t_span=time_span,

问题二

建立问题二的模型需要考虑直升机在飞行中的滚转、俯仰和偏航时的力矩模型。我们可以采用欧拉角（Euler angles）来描述直升机的姿态，其中包括滚转角（roll angle）、俯仰角（pitch angle）和偏航角（yaw angle）。

1. 动力学方程

在飞行器的滚转、俯仰和偏航运动中，牛顿-欧拉方程可以表示为：

$I_x\cdot \dot{p}=(I_y-I_z)\cdot q\cdot r+{\tau }_{\text{control-roll}}$

$I_y\cdot \dot{q}=(I_z-I_x)\cdot p\cdot r+{\tau }_{\text{control-pitch}}$

$I_z\cdot \dot{r}=(I_x-I_y)\cdot p\cdot q+{\tau }_{\text{control-yaw}}$

其中：

• $I_x,I_y,I_z$分别是绕X、Y、Z轴的惯性矩；
• $p,q,r$是滚转、俯仰和偏航角速度；
• ${\tau }_{\text{control-roll}},{\tau }_{\text{control-pitch}},{\tau }_{\text{control-yaw}}$是分别由滚转、俯仰和偏航控制产生的力矩。

2. 各组件的力矩

考虑直升机上的各个组件对力矩的影响：

滚转控制：

${\tau }_{\text{control-roll}}=k_{\text{elevator}}\cdot {\delta }_{\text{elevator}}$

俯仰控制：

${\tau }_{\text{control-pitch}}=k_{\text{elevator}}\cdot {\delta }_{\text{elevator}}$

偏航控制：

${\tau }_{\text{control-yaw}}=k_{\text{rudder}}\cdot {\delta }_{\text{rudder}}$

其中：

• $k_{\text{elevator}},k_{\text{rudder}}$是与升降舵和方向舵的控制效率有关的常数；
• ${\delta }_{\text{elevator}},{\delta }_{\text{rudder}}$是对应的控制输入。

3. 数学模型

将以上方程整合，得到关于滚转、俯仰和偏航的动力学方程：

$I_x\cdot \dot{p}=(I_y-I_z)\cdot q\cdot r+k_{\text{elevator}}\cdot {\delta }_{\text{elevator}}$

$I_y\cdot \dot{q}=(I_z-I_x)\cdot p\cdot r+k_{\text{elevator}}\cdot {\delta }_{\text{elevator}}$

KaTeX parse error: Expected '}', got 'EOF' at end of input: …_{\text{rudder}

4. 数值模拟

使用数值模拟工具，例如ODE求解器，模拟直升机在给定初始条件和控制输入的情况下，滚转、俯仰和偏航角随时间的演化。

5. 代码

import numpy as np
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt

# 定义直升机姿态角动力学方程
def attitude_dynamics(t, y, elevator_input, rudder_input):
    # y: [roll, pitch, yaw, roll_rate, pitch_rate, yaw_rate]
    roll, pitch, yaw, roll_rate, pitch_rate, yaw_rate = y
    
    # 参数设定
    I_x, I_y, I_z = 1000.0, 1500.0, 2000.0  # 代表绕各轴的惯性矩，实际应根据直升机参数调整
    k_elevator, k_rudder = 2.0, 2.0  # 控制效率常数
    
    # 控制输入
    delta_elevator = elevator_input(t)
    delta_rudder = rudder_input(t)
    
    # 力矩模型
    control_roll_torque = k_elevator * delta_elevator
    control_pitch_torque = k_elevator * delta_elevator
    control_yaw_torque = k_rudder * delta_rudder
    
    # 动力学方程
    roll_acceleration = ((I_y - I_z) * pitch_rate * yaw_rate + control_roll_torque) / I_x
    pitch_acceleration = ((I_z - I_x) * roll_rate * yaw_rate + control_pitch_torque) / I_y
    yaw_acceleration = ((I_x - I_y) * roll_rate * pitch_rate + control_yaw_torque) / I_z

    return [roll_rate, pitch_rate, yaw_rate, roll_acceleration, pitch_acceleration, yaw_acceleration]

# 初始条件
initial_conditions = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0]  # 初始角度为0，初始角速度为0

# 时间范围
time_span = (0, 20)

# 控制输入函数（示例，可以根据需要更改）
def elevator_input(t):
#见完整代码

问题三

问题三要求基于直升机在低速和高速飞行中的机动特性，设计每个组件的机动幅度，以满足直升机在水平飞行任务下的需求。这里我们可以采用最优控制理论，具体分为两个阶段：低速飞行和高速飞行。

1. 低速飞行

a. 满足水平飞行任务

在低速飞行时，水平飞行任务是维持零姿态角。我们可以使用最优控制理论，通过调整每个组件的控制输入，使得直升机在低速水平飞行时保持零姿态角。

b. 设计每个组件的机动幅度

通过调整每个组件的控制输入，例如调整升降舵和方向舵的偏转角，以实现水平飞行任务。这涉及到设计合适的控制输入幅度和频率。

2. 高速飞行

a. 满足水平飞行任务

在高速飞行时，同样需要满足水平飞行任务，即保持零姿态角。

b. 设计每个组件的机动幅度

在高速飞行时，由于机动要求可能更为复杂，我们需要调整每个组件的控制输入，以满足水平飞行任务。这可能涉及到更大的控制输入幅度和更高的频率。

3. 数学建模

使用最优控制理论，可以建立一个优化问题，目标是最小化控制输入的幅度，同时满足水平飞行任务的约束。这可以通过设定一个性能指标，例如最小化控制输入的平方和。

4. 代码实现

具体实现涉及使用最优控制理论的数学工具和算法，例如 Pontryagin最大值原理等。这可能需要使用专门的最优化库，如SciPy中的优化工具。

import numpy as np
from scipy.integrate import solve_ivp
import scipy.optimize as optimize
import matplotlib.pyplot as plt

# 定义目标函数
def cost_function(inputs):
    # inputs: 控制输入向量
    # 模拟直升机在给定控制输入下的姿态角演化
    # 这里假设存在一个模拟函数 simulate_attitude(inputs)，返回姿态角随时间的演化
    simulated_attitude = simulate_attitude(inputs)
    
    # 计算性能指标，例如姿态角偏差的平方和
    target_attitude = np.zeros_like(simulated_attitude)  # 零姿态角作为目标
    cost = np.sum((simulated_attitude - target_attitude)**2)
    
    return cost

# 定义约束条件，例如水平飞行任务的约束条件
def horizontal_flight_constraint(inputs):
    # inputs: 控制输入向量
    # 返回水平飞行任务的约束条件，可以是姿态角的某个范围或其他条件
    # 这里简化为控制输入的幅度不能超过某个阈值
    return np.sum(np.abs(inputs)) - max_input_threshold

# 模拟函数，用于计算在给定控制输入下直升机的姿态角演化
def simulate_attitude(inputs):
    # 使用ODE求解器模拟直升机姿态角的演化
    # 这里假设存在一个 attitude_dynamics 函数表示直升机的姿态动力学方程
    solution = solve_ivp(attitude_dynamics, time_span, initial_conditions, args=(inputs,), method='RK45', dense_output=True)
    
    # 获取模拟结果
    simulated_attitude = solution.sol(t_eval)
    
    return simulated_attitude

# 初始控制输入猜测值
initial_guess = np.array([0.1, 0.1, 0.1, ...])

# 控制输入的约束条件
max_input_threshold = 1.0
constraints = {'type': 'eq', 'fun': horizontal_flight_constraint}

# 设置优化问题
optimization_problem = optimize.minimize(cost_function, initial_guess, constraints=constraints)

# 获取优化结果
optimal_inputs = optimization_problem.x

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

1. 目标：

设计每个控制输入的动态值，使得直升机在20秒内均匀加速，飞行速度从80 m/s增加到180 m/s，同时保持水平飞行（零姿态角）。

2. 建模思路：

a. 数学建模：

• 使用最优控制理论，考虑直升机姿态角随时间的演化。
• 定义目标函数，例如姿态角偏差的平方和，以表示性能指标。
• 建立姿态动力学方程，考虑控制输入对姿态角的影响。

b. 控制输入动态调整：

• 动态调整每个控制输入，使直升机实现前向加速度。
• 根据加速度需求，调整控制输入的幅度和频率。

3. 详细思路：

1. 建立目标函数：

   • 定义目标函数，例如$J={\sum }_{i=1}^N({\theta }_i-{\theta }_{i,\text{target}}{)}^2$，其中 ${\theta }_i$ 表示姿态角，$N$ 表示时间步数，${\theta }_{i,\text{target}}$表示目标姿态角。

2. 模拟姿态角演化：

   • 使用ODE求解器模拟直升机姿态角的演化。
   • 通过姿态动力学方程，考虑控制输入对姿态角的影响。

3. 控制输入动态调整：

   • 在模拟过程中，动态调整每个控制输入，以满足前向加速度的需求。
   • 根据加速度目标，动态地调整控制输入的幅度和频率。

4. 优化问题求解：

   • 利用最优化算法，例如SciPy库中的优化工具，求解目标函数的最小值。
   • 初始控制输入猜测值可以通过启发式方法或者根据问题背景设定。

5. 模拟最优控制输入下的姿态角演化：

   • 使用求解得到的最优控制输入，模拟直升机在加速机动任务中的姿态角演化。

6. 结果可视化：

   • 绘制模拟结果，展示直升机姿态角随时间的变化。
   • 检查结果是否满足前向加速度和水平飞行任务。

4. 代码实现：

import numpy as np
from scipy.integrate import solve_ivp
import scipy.optimize as optimize
import matplotlib.pyplot as plt

# 假设直升机的动力学方程和控制特性如下

# 控制输入：每个控制输入表示直升机不同部分的控制量
# 例如：inputs = [main_rotor_pitch, tail_rotor_pitch, ...]

# 目标函数
def cost_function(inputs):
    # 模拟直升机在给定控制输入下的姿态角演化
    simulated_attitude = simulate_attitude(inputs)
    
    # 计算性能指标，例如姿态角偏差的平方和
    target_attitude = np.zeros_like(simulated_attitude)  # 零姿态角作为目标
    cost = np.sum((simulated_attitude - target_attitude)**2)
    
    return cost

# 模拟函数
def simulate_attitude(inputs):
    # 使用ODE求解器模拟直升机姿态角的演化
    solution = solve_ivp(attitude_dynamics, time_span, initial_conditions, args=(inputs,), method='RK45', dense_output=True)
    
    # 获取模拟结果
    simulated_attitude = solution.sol(t_eval)
    
    return simulated_attitude

# 姿态动力学方程
def attitude_dynamics(t, y, inputs):
    # y: [roll, pitch, yaw, roll_rate, pitch_rate, yaw_rate]
    # 根据控制输入计算每个控制力和力矩
    
    # 控制输入的动态调整
    # 这里可以根据加速度的需求动态调整控制输入，实现前向加速度
    
    # 姿态动力学方程
    # 这里假设存在合适的动力学方程，需要根据直升机的具体参数和特性进行实现
    roll_acceleration = 0.0  # 根据实际情况计算
    pitch_acceleration = 0.0  # 根据实际情况计算
    yaw_acceleration = 0.0  # 根据实际情况计算
    
    return [roll_rate, pitch_rate, yaw_rate, roll_acceleration, pitch_acceleration, yaw_acceleration]

# 初始控制输入猜测值
initial_guess = np.array([0.1, 0.1, 0.1, ...])

# 时间范围
time_span = (0, 20)

# 优化问题求解
optimization_problem = optimize.minimize(cost_function, initial_guess, method='SLSQP')

# 获取优化结果
optimal_inputs = optimization_problem.x

# 模拟最优控制输入下的姿态角演化
simulated_attitude = simulate_attitude(optimal_inputs)

# 绘制结果
plt.figure(figsize=(12, 8))
plt.plot(t_eval, np.degrees(simulated_attitude[0]), label='Roll Angle (deg)')
plt.plot(t_eval, np.degrees(simulated_attitude[1]), label='Pitch Angle (deg)')
plt.plot(t_eval, np.degrees(simulated_attitude[2]), label='Yaw Angle (deg)')

更多完整的思路和代码看这里哦：
2023 第九届数维杯A题 建模解析，小鹿学长带队