1. 从“空中芭蕾”到数学建模：一个动力学问题的拆解

看到“空中芭蕾”这个题目，很多同学的第一反应可能是艺术体操或者花样滑冰。但在2025年数维杯数学建模的赛场上，它指向的是一个典型的、充满挑战的动力学与控制问题。这类问题往往描述一个或多个物体（比如飞行器、机器人、甚至是一个复杂的多连杆系统）在三维空间中的运动，要求其轨迹不仅满足物理规律，还要具备某种“优美”或“高效”的特性——就像芭蕾舞演员的动作一样，既精准又富有表现力。

我参加过不少数学建模竞赛，也带过很多队伍。我发现，面对A题这种开放性较强的题目，最大的障碍不是数学或编程，而是如何将一段充满想象力的文字描述，转化成一个结构清晰、可量化、可求解的数学模型。“空中芭蕾”听起来很抽象，但它本质上逃不开几个核心要素：研究对象（谁在跳舞？）、约束条件（舞台有多大？规则是什么？）、目标函数（怎样才算跳得好？）、以及环境与干扰（有没有风？地面滑不滑？）。这篇内容，我就结合这个题目，把从审题到代码落地的完整思考过程和实操细节拆开揉碎了讲给你听。无论你是第一次参赛的新手，还是想提升解题层次的老手，希望这些从实战中踩坑得来的经验，能帮你少走弯路。

我们的目标很明确：第一，彻底读懂“空中芭蕾”背后隐藏的物理与数学问题；第二，建立一套有效的建模框架；第三，用Python将模型“实现”出来，得到可视化结果与定量分析；第四，也是最重要的，将这些过程整理成逻辑严谨、可复现的论文。你会发现，编程（Python）只是工具，真正的灵魂在于你如何定义问题、做出假设、并设计求解策略。

2. 问题重述与核心概念界定：把“舞蹈”翻译成“方程”

拿到赛题，千万别急着翻算法书或者写代码。第一步，也是决定成败的一步，是精确地重述问题。我们以“空中芭蕾”为假想题，来模拟这个过程。

2.1 拆解题目关键词

假设原题描述大致为：“设计一个控制系统，使一个多旋翼无人机（或一个简化质点）在三维空间中完成一段复杂的特技飞行轨迹。该轨迹需满足：1. 整体时间T内完成；2. 途经若干个指定的关键点（位置、姿态）；3. 飞行平滑，加速度变化连续；4. 能量消耗尽可能小。试建立数学模型，规划最优轨迹，并评估其性能。”

从这个假想描述中，我们可以剥离出以下核心建模要素：

• 系统模型：我们控制的对象是什么？它的动力学方程是怎样的？对于多旋翼无人机，这是一个典型的欠驱动系统（只有4个控制输入：总升力和三个姿态角扭矩，却要控制6个自由度：位置和姿态）。我们可能需要对其进行简化，例如在轨迹规划层，先忽略复杂的姿态动力学，将其视为一个可控质点，其加速度作为直接控制量。这在建模中称为“微分平坦”特性利用，是简化问题的关键。
• 状态变量：描述系统“状态”的最小变量集。对于质点模型，状态就是位置(x, y, z)和速度(v_x, v_y, v_z)。对于刚体模型，还需加入姿态（如四元数[qw, qx, qy, qz]）和角速度。
• 控制变量：我们可以直接操纵的量。对于质点模型，可以是加速度(a_x, a_y, a_z)；对于真实无人机，则是四个电机的转速。
• 路径约束：飞行中必须时刻满足的条件。例如，速度不能超过v_max，加速度（反映电机推力极限）不能超过a_max，位置必须处于安全空域内。
• 边界条件：起始点和终止点的状态。通常给定起始位置速度、终止位置速度（例如，都从静止开始并结束于静止）。
• 中间点约束：轨迹必须精确穿过或以一定精度接近的某些关键点（“舞步”的节点）。可能还要求在到达这些点时具有特定的速度方向或大小。
• 目标函数（性能指标）：用来评价轨迹“好坏”的数学标准。“能量消耗尽可能小”是一个常见指标。对于无人机，能量消耗通常与控制量的平方（反映电机推力）的积分成正比，即最小化∫ (u_x² + u_y² + u_z²) dt。也可能追求时间最短min T，或者追求平滑度（控制量的变化率，即加加速度jerk最小）。

2.2 做出合理假设

数学建模离不开合理的假设，这是将复杂现实问题变为可解数学问题的桥梁。针对“空中芭蕾”，我们可能会做如下假设：

1. 假设空气阻力与速度成正比：这是一个常见的线性简化，比二次型阻力更容易处理。公式可写为F_drag = -k * v。如果题目强调高速或精确，则需考虑二次型阻力。
2. 假设无人机为质点：专注于轨迹规划，暂不考虑姿态动力学。这样控制变量直接是加速度，问题简化为一个带约束的质点运动规划。
3. 假设重力加速度恒定：g = 9.8 m/s²，方向垂直向下。
4. 假设关键点之间无碰撞风险：即只规划单机轨迹，或多机间有预设的安全间隔。 这些假设必须在论文中明确列出，并简要说明其合理性。它们定义了模型的适用范围。

2.3 确定建模与求解思路

明确了问题和假设，接下来要选择“打法”。对于轨迹优化问题，主流思路有两类：

• 最优控制理论：将问题表述为庞特里亚金最小值原理（PMP）下的两点边值问题，或者直接使用哈密顿-雅可比-贝尔曼方程（HJB）。这种方法理论优美，能得到全局最优的必要条件，但对于复杂约束和非线性系统，解析求解极其困难，数值求解（如打靶法）对初值敏感。
• 数值优化方法：将连续时间问题离散化，转化为一个大规模的、有约束的数值优化问题。这是目前实践中最主流、最强大的方法。具体来说，就是把总时间T分成N个小段，每一段的状态和控制量都变成优化变量，动力学方程转化为约束条件，然后用非线性规划求解器求解。

我们将采用数值优化方法，因为它更通用，能方便地处理各种复杂约束，并且有成熟的Python工具包（如CasADi）支持。思路框架如下：

1. 离散化：将时间[0, T]等分为N个区间，得到N+1个时间节点。
2. 定义决策变量：每个时间节点上的状态变量（位置、速度）和控制变量（加速度）都成为优化变量。总共约有(6 + 3) * (N+1)个标量变量（质点模型）。
3. 构造约束：
   • 动力学约束：使用数值积分方法（如欧拉法、中点法、龙格-库塔法）将连续动力学方程dv/dt = a, dx/dt = v转化为相邻离散状态变量之间的等式约束。例如，使用欧拉法：x[k+1] = x[k] + v[k] * dt,v[k+1] = v[k] + a[k] * dt。
   • 路径约束：v_min <= |v[k]| <= v_max,|a[k]| <= a_max。
   • 边界约束：x[0], v[0], x[N], v[N]等于给定的值。
   • 中间点约束：对于指定的关键时间点k_j，要求x[k_j] = x_target_j。
4. 构造目标函数：例如，最小化控制力平方和∑ |a[k]|² * dt。
5. 调用求解器：将上述问题输入非线性规划求解器（如IPOPT）进行求解。

3. 基于CasADi的轨迹优化Python实现详解

理论清晰后，我们来动手实现。Python中，CasADi是一个专门用于优化和最优控制的强大框架，它支持自动微分，并能高效地接口IPOPT等求解器。下面，我们一步步构建一个完整的“空中芭蕾”轨迹优化程序。

3.1 环境准备与问题参数定义

首先，确保你的Python环境安装了必要的库。在终端或命令提示符中执行：

pip install casadi numpy matplotlib

接下来，开始编写代码。我们首先定义问题的物理参数和优化参数。

import casadi as ca import numpy as np import matplotlib.pyplot as plt from mpl_toolkits.mplot3d import Axes3D # ==================== 1. 问题参数定义 ==================== # 物理参数 g = 9.8 # 重力加速度，m/s^2 mass = 1.0 # 无人机质量，kg k_drag = 0.1 # 线性阻尼系数，kg/s (假设) # 轨迹参数 T = 10.0 # 总飞行时间，秒 N = 100 # 离散化网格数量，即将T分成N段，有N+1个点 dt = T / N # 时间间隔 # 状态与控制的边界约束 v_max = 5.0 # 最大速度，m/s a_max = 3.0 # 最大加速度，m/s^2 # 起始和终止状态 x_start = [0, 0, 0] # 起始位置 (x, y, z) v_start = [0, 0, 0] # 起始速度 x_goal = [10, 5, 3] # 目标位置 v_goal = [0, 0, 0] # 目标速度 # 定义中间关键点（“舞步”节点） # 格式：[时间比例 (0~1), x, y, z] waypoints = [ [0.2, 2, 1, 2], [0.5, 5, 3, 4], [0.8, 8, 4, 2.5] ]

注意：这里我们使用线性阻尼模型，并将关键点定义为时间比例。在实际比赛中，题目可能要求精确的时间点或顺序经过某些点，需要相应调整约束形式。

3.2 构建优化问题：变量、约束与目标

这是最核心的部分。我们使用CasADi的符号变量来构建优化问题。

# ==================== 2. 建立优化问题 ==================== # 定义优化变量容器 opti = ca.Opti() # 决策变量：每个时间步的状态和控制量 # 状态：位置 (3维) + 速度 (3维) = 6维 X = opti.variable(6, N+1) # 维度： (6状态) x (N+1个时间点) # 控制：加速度 (3维) U = opti.variable(3, N) # 维度： (3控制) x (N个控制区间) # 将状态变量拆分为位置和速度，便于书写 pos = X[0:3, :] # 3 x (N+1) vel = X[3:6, :] # 3 x (N+1) acc = U # 3 x N # -------------------- 2.1 动力学约束（离散化） -------------------- # 使用简单的欧拉前向差分进行离散化 for k in range(N): # 位置更新: pos[k+1] = pos[k] + vel[k] * dt opti.subject_to(pos[:, k+1] == pos[:, k] + vel[:, k] * dt) # 速度更新: vel[k+1] = vel[k] + (acc[:, k] - k_drag*vel[:, k]/mass - [0,0,g]) * dt # 注意：这里包含了阻尼力和重力。控制加速度acc是总加速度的一部分。 # 实际动力学: m * dv/dt = F_control - k_drag * v - m*g*e_z # 所以 dv/dt = (F_control/m) - (k_drag/m)*v - g*e_z # 令我们的控制变量 U (acc) 代表 (F_control/m)，即单位质量的推力。 drag_force = (k_drag / mass) * vel[:, k] gravity = ca.vertcat(0, 0, g) opti.subject_to(vel[:, k+1] == vel[:, k] + (acc[:, k] - drag_force - gravity) * dt) # -------------------- 2.2 边界条件约束 -------------------- # 初始状态 opti.subject_to(pos[:, 0] == x_start) opti.subject_to(vel[:, 0] == v_start) # 终端状态 opti.subject_to(pos[:, N] == x_goal) opti.subject_to(vel[:, N] == v_goal) # -------------------- 2.3 路径约束 -------------------- # 速度幅值约束 for k in range(N+1): v_norm_squared = vel[0, k]**2 + vel[1, k]**2 + vel[2, k]**2 opti.subject_to(v_norm_squared <= v_max**2) # 控制加速度幅值约束 for k in range(N): a_norm_squared = acc[0, k]**2 + acc[1, k]**2 + acc[2, k]**2 opti.subject_to(a_norm_squared <= a_max**2) # -------------------- 2.4 中间点（关键点）约束 -------------------- for wp in waypoints: t_ratio, x_t, y_t, z_t = wp k_idx = int(round(t_ratio * N)) # 找到对应的时间节点索引 # 确保索引在有效范围内 k_idx = max(0, min(N, k_idx)) opti.subject_to(pos[0, k_idx] == x_t) opti.subject_to(pos[1, k_idx] == y_t) opti.subject_to(pos[2, k_idx] == z_t) # 如果需要，也可以在这里添加速度方向约束，例如要求到达该点时速度水平： # opti.subject_to(vel[2, k_idx] == 0) # -------------------- 2.5 目标函数 -------------------- # 最小化控制努力（加速度平方和），这是一个常见的能量消耗代理指标 objective = 0 for k in range(N): objective += ca.sumsqr(acc[:, k]) * dt # 积分近似为求和 opti.minimize(objective) # -------------------- 2.6 提供初始猜测 -------------------- # 一个好的初始猜测能显著帮助求解器收敛。这里我们用简单的线性插值。 # 位置初始猜测：从起点到终点的直线 for i in range(3): opti.set_initial(pos[i, :], np.linspace([x_start[i], x_goal[i]], N+1)) # 速度初始猜测：设为0 opti.set_initial(vel, 0) # 控制初始猜测：设为0 opti.set_initial(acc, 0)

关键点解析：

1. 动力学离散化：这里使用了显式欧拉法，简单但精度一般。对于更高精度的要求，可以使用中点法 (RK2) 或RK4。CasADi也内置了integrator函数来处理更复杂的连续动力学。
2. 约束添加：opti.subject_to()是添加约束的核心函数。注意我们对每个时间点循环添加约束，这可能会产生大量约束条件，但CasADi和IPOPT能高效处理稀疏结构。
3. 目标函数：最小化控制量的平方和（即∫ ||a(t)||² dt），这通常对应于最小化能量消耗或控制力波动，能使轨迹更平滑。
4. 初始猜测：对于非线性优化问题，初始值至关重要。线性插值是一个简单有效的选择。如果问题非凸（例如存在障碍物），可能需要更复杂的初始化策略，甚至使用多起点优化。

3.3 求解与结果提取

问题构建完成后，就可以调用求解器了。

# ==================== 3. 求解优化问题 ==================== # 选择求解器（IPOPT是CasADi默认集成的强大开源求解器） opti.solver('ipopt') # 尝试求解 try: sol = opti.solve() print("优化求解成功！") except RuntimeError as e: print(f"求解失败: {e}") # 如果失败，可以尝试查看不可行的约束 # 或者使用调试模式：opti.debug.show_infeasibilities() # 这里我们退出或进行降级处理 exit(1) # 提取最优解 pos_opt = sol.value(pos) vel_opt = sol.value(vel) acc_opt = sol.value(acc) time_axis = np.linspace(0, T, N+1)

3.4 可视化与结果分析

将优化得到的轨迹画出来，是验证模型和呈现结果的关键。

# ==================== 4. 结果可视化 ==================== fig = plt.figure(figsize=(18, 10)) # 1. 3D轨迹图 ax1 = fig.add_subplot(231, projection='3d') ax1.plot(pos_opt[0, :], pos_opt[1, :], pos_opt[2, :], 'b-', linewidth=2, label='Optimal Trajectory') ax1.scatter(pos_opt[0, 0], pos_opt[1, 0], pos_opt[2, 0], c='g', s=100, marker='o', label='Start') ax1.scatter(pos_opt[0, -1], pos_opt[1, -1], pos_opt[2, -1], c='r', s=100, marker='^', label='Goal') for wp in waypoints: _, x, y, z = wp ax1.scatter(x, y, z, c='m', s=80, marker='s', label='Waypoint' if wp == waypoints[0] else "") ax1.set_xlabel('X [m]') ax1.set_ylabel('Y [m]') ax1.set_zlabel('Z [m]') ax1.set_title('3D Optimal Trajectory') ax1.legend() ax1.grid(True) # 2. 位置-时间曲线 ax2 = fig.add_subplot(234) ax2.plot(time_axis, pos_opt[0, :], label='X') ax2.plot(time_axis, pos_opt[1, :], label='Y') ax2.plot(time_axis, pos_opt[2, :], label='Z') for wp in waypoints: t_ratio, x, y, z = wp t = t_ratio * T ax2.axvline(x=t, color='gray', linestyle='--', alpha=0.5) ax2.set_xlabel('Time [s]') ax2.set_ylabel('Position [m]') ax2.set_title('Position vs. Time') ax2.legend() ax2.grid(True) # 3. 速度-时间曲线 ax3 = fig.add_subplot(235) vel_norm = np.linalg.norm(vel_opt, axis=0) ax3.plot(time_axis, vel_opt[0, :], label='Vx') ax3.plot(time_axis, vel_opt[1, :], label='Vy') ax3.plot(time_axis, vel_opt[2, :], label='Vz') ax3.plot(time_axis, vel_norm, 'k--', label='Norm', linewidth=2) ax3.axhline(y=v_max, color='r', linestyle='--', label=f'V_max ({v_max} m/s)') ax3.set_xlabel('Time [s]') ax3.set_ylabel('Velocity [m/s]') ax3.set_title('Velocity vs. Time') ax3.legend() ax3.grid(True) # 4. 控制加速度-时间曲线 ax4 = fig.add_subplot(236) acc_norm = np.linalg.norm(acc_opt, axis=0) # 注意控制量只有N个点，时间轴需要调整 time_control = np.linspace(0, T, N+1)[:-1] + dt/2 # 取区间中点代表该区间控制量 ax4.plot(time_control, acc_opt[0, :], label='Ax') ax4.plot(time_control, acc_opt[1, :], label='Ay') ax4.plot(time_control, acc_opt[2, :], label='Az') ax4.plot(time_control, acc_norm, 'k--', label='Norm', linewidth=2) ax4.axhline(y=a_max, color='r', linestyle='--', label=f'A_max ({a_max} m/s²)') ax4.set_xlabel('Time [s]') ax4.set_ylabel('Control Acceleration [m/s²]') ax4.set_title('Control Input vs. Time') ax4.legend() ax4.grid(True) plt.tight_layout() plt.show() # 打印一些关键性能指标 print("\n=== 轨迹性能指标 ===") print(f"总控制努力 (∫||a||² dt): {sol.value(objective):.4f}") print(f"最大速度: {vel_norm.max():.4f} m/s") print(f"最大控制加速度: {acc_norm.max():.4f} m/s²") print(f"是否满足速度约束: {vel_norm.max() <= v_max + 1e-6}") print(f"是否满足加速度约束: {acc_norm.max() <= a_max + 1e-6}")

运行这段代码，你将得到四张图：3D轨迹图、位置-时间曲线、速度-时间曲线和控制量-时间曲线。从图中可以直观地看到，无人机是如何平滑地飞过所有关键点，并且速度和加速度都严格限制在约束范围内。控制量曲线通常比较“饱满”，这意味着求解器在充分利用系统的能力来优化目标。

4. 模型扩展、论文写作要点与实战避坑指南

一个基础的模型跑通后，真正的挑战在于如何让它更贴近实际、更鲁棒，以及如何将整个工作清晰地呈现在论文中。

4.1 模型进阶与扩展方向

上述模型是一个起点。针对“空中芭蕾”或类似赛题，你可以从以下几个方向深化模型，这往往是论文的加分项：

• 更精确的动力学模型：将质点模型升级为刚体模型。引入姿态动力学（四元数或欧拉角表示），控制量变为电机的推力。这会显著增加问题的非线性程度和变量规模，但能研究姿态与轨迹的耦合效应。
• 考虑障碍物与避碰：在空域中添加圆柱体、长方体等障碍物。这需要在路径约束中添加非凸约束（例如，位置到障碍物表面的距离大于安全阈值d_min）。处理非凸约束是难点，可能需要引入整数变量（混合整数规划）或使用序列凸优化（SCP）等高级技巧。
• 多智能体协同“芭蕾”：规划多架无人机的轨迹，要求它们保持队形、避免碰撞。这会引入大量的相对位置约束，问题规模呈组合爆炸。一种策略是先规划质心轨迹，再规划相对运动。
• 不确定性处理：考虑风扰、模型参数不确定性。这可以引入鲁棒优化或随机优化的框架，目标函数变为最小化期望成本或最坏情况成本。
• 微分平坦性应用：对于许多机器人系统（如四旋翼），其轨迹规划可以在输出空间（位置、偏航角）进行，然后通过微分平坦变换反推状态和控制量。这能极大简化问题，将轨迹参数化为多项式或B样条曲线，只优化曲线参数。

4.2 数学建模论文写作的核心结构

论文是将你的思想呈现给评委的唯一载体。结构清晰、逻辑自洽至关重要。一篇完整的数模论文通常包含：

1. 摘要：重中之重！用300-500字概括整个工作。必须包含：问题重述、你的建模思路、所用方法、主要结果、结论与特色。避免细节，突出整体逻辑和创新点。最后写，但最先被阅读。
2. 问题重述与分析：用自己的语言精炼地复述问题，并进行分析，指出问题的关键点、难点和可能的解决路径。展示你对问题的深刻理解。
3. 模型假设与符号说明：清晰列出所有假设，并说明其合理性。制作一个符号表，列出所有主要变量、符号及其含义和单位。
4. 模型建立与求解：这是论文的躯干。
   • 模型准备：阐述建模框架（如最优控制、数值优化）。
   • 模型建立：详细推导状态方程、约束条件和目标函数。公式要编号，解释要清晰。
   • 模型求解：说明你采用的求解算法（如直接转录法+IPOPT）、离散化方法、以及如何处理特定约束（如中间点约束）。可以配上算法流程图。
   • 参数设置：给出模型中所有参数的值及其来源（题目给定、合理假设或参考文献）。
5. 模型求解与结果分析：
   • 仿真环境：说明你的编程语言、工具包和计算平台。
   • 基准案例：展示一个标准参数下的完整结果（就像我们上面代码跑出来的那样），配上精心设计的图表和文字分析。分析要指出：“由图X可见，轨迹平滑穿过所有关键点...”、“表Y显示，最大速度未超过约束...”、“这表明模型能有效生成可行轨迹...”。
   • 灵敏度分析：这是体现模型稳健性和你分析能力的关键！改变关键参数（如v_max,a_max, 关键点位置，总时间T），观察目标函数和轨迹形态如何变化。分析其物理意义和工程启示。
   • 模型对比/扩展：如果时间允许，将你的模型与一个简单模型（如不考虑能量的最短路径）对比，或者展示扩展模型（如加入障碍物）的结果。
6. 模型评价与推广：客观评价模型的优点（如考虑全面、求解高效、结果合理）和缺点（如假设简化、未考虑XX因素）。提出模型的改进方向和在更广领域的应用前景。
7. 参考文献：规范引用。
8. 附录：可以放置核心代码（不必全部，关键片段即可）。

4.3 实战中的常见“坑”与应对策略

结合我自己的经验，分享几个最容易出问题的地方：

• 求解器不收敛或求解速度慢：
  • 原因：初始猜测太差、约束相互冲突、问题尺度或数值条件太差。
  • 对策：1) 提供更好的初始值（如线性插值）。2) 逐步增加问题复杂度，先解一个简化问题（如去掉中间点约束），用其解作为完整问题的初始猜测。3) 检查约束是否可能矛盾（例如，给定的总时间T太短，无法在速度限制下到达终点）。4) 缩放变量，让所有决策变量和约束的量级都在1附近（例如，位置除以10，速度除以5），这能极大改善求解器的数值稳定性。
• 结果违反约束：
  • 原因：IPOPT等求解器允许微小的约束违反（容差tol）。有时是离散化误差或数值误差导致。
  • 对策：1) 检查违反的程度，如果远小于容差（如1e-6），可以认为是数值噪声，在论文中说明即可。2) 如果违反较大，尝试收紧求解器容差（opti.solver('ipopt', {'tol':1e-8})），或增加离散化网格数量N。
• 代码调试困难：
  • 对策：善用opti.debug。当求解失败时，opti.debug.show_infeasibilities()可以高亮最可能违反的约束，是定位问题的神器。另外，将大问题分解为小问题分段测试。
• 论文图表不专业：
  • 对策：图表是论文的门面。确保所有图表都有清晰的标题、坐标轴标签（带单位）、图例。线型、颜色要区分明显。3D图要选择合适的视角，能清晰展示轨迹。避免使用默认的过于花哨的样式，保持简洁、学术化。Matplotlib的plt.style.use('seaborn-v0_8-whitegrid')是一个不错的起点。
• 时间管理失控：
  • 对策：三天比赛，第一天必须完成问题分析、初步建模和基础代码框架。第二天上午要得到第一个可运行的结果，下午进行深入分析和模型改进/扩展。第三天全天用于论文写作、润色和制作图表。编程和建模要并行，不要等模型完美了再写代码，也不要代码跑通了再开始写论文。论文写作是持续的过程。

最后，我想强调的是，“空中芭蕾”这类题目考察的远不止编程和数学。它考察的是你将模糊需求转化为精确模型的能力、在复杂约束下寻找解决方案的创造力，以及用严谨文字和图表展示工作的沟通力。上面的代码和框架提供了一个坚实的起点，但真正的亮点来自于你对问题的独特见解和深入分析。多思考“如果...会怎样？”，并尝试用你的模型去回答，你的论文就会脱颖而出。