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 做出合理假设数学建模离不开合理的假设这是将复杂现实问题变为可解数学问题的桥梁。针对“空中芭蕾”我们可能会做如下假设假设空气阻力与速度成正比这是一个常见的线性简化比二次型阻力更容易处理。公式可写为F_drag -k * v。如果题目强调高速或精确则需考虑二次型阻力。假设无人机为质点专注于轨迹规划暂不考虑姿态动力学。这样控制变量直接是加速度问题简化为一个带约束的质点运动规划。假设重力加速度恒定g 9.8 m/s²方向垂直向下。假设关键点之间无碰撞风险即只规划单机轨迹或多机间有预设的安全间隔。 这些假设必须在论文中明确列出并简要说明其合理性。它们定义了模型的适用范围。2.3 确定建模与求解思路明确了问题和假设接下来要选择“打法”。对于轨迹优化问题主流思路有两类最优控制理论将问题表述为庞特里亚金最小值原理PMP下的两点边值问题或者直接使用哈密顿-雅可比-贝尔曼方程HJB。这种方法理论优美能得到全局最优的必要条件但对于复杂约束和非线性系统解析求解极其困难数值求解如打靶法对初值敏感。数值优化方法将连续时间问题离散化转化为一个大规模的、有约束的数值优化问题。这是目前实践中最主流、最强大的方法。具体来说就是把总时间T分成N个小段每一段的状态和控制量都变成优化变量动力学方程转化为约束条件然后用非线性规划求解器求解。我们将采用数值优化方法因为它更通用能方便地处理各种复杂约束并且有成熟的Python工具包如CasADi支持。思路框架如下离散化将时间[0, T]等分为N个区间得到N1个时间节点。定义决策变量每个时间节点上的状态变量位置、速度和控制变量加速度都成为优化变量。总共约有(6 3) * (N1)个标量变量质点模型。构造约束动力学约束使用数值积分方法如欧拉法、中点法、龙格-库塔法将连续动力学方程dv/dt a, dx/dt v转化为相邻离散状态变量之间的等式约束。例如使用欧拉法x[k1] x[k] v[k] * dt,v[k1] 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。构造目标函数例如最小化控制力平方和∑ |a[k]|² * dt。调用求解器将上述问题输入非线性规划求解器如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段有N1个点 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, N1) # 维度 (6状态) x (N1个时间点) # 控制加速度 (3维) U opti.variable(3, N) # 维度 (3控制) x (N个控制区间) # 将状态变量拆分为位置和速度便于书写 pos X[0:3, :] # 3 x (N1) vel X[3:6, :] # 3 x (N1) acc U # 3 x N # -------------------- 2.1 动力学约束离散化 -------------------- # 使用简单的欧拉前向差分进行离散化 for k in range(N): # 位置更新: pos[k1] pos[k] vel[k] * dt opti.subject_to(pos[:, k1] pos[:, k] vel[:, k] * dt) # 速度更新: vel[k1] 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[:, k1] 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(N1): 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]], N1)) # 速度初始猜测设为0 opti.set_initial(vel, 0) # 控制初始猜测设为0 opti.set_initial(acc, 0)关键点解析动力学离散化这里使用了显式欧拉法简单但精度一般。对于更高精度的要求可以使用中点法 (RK2) 或RK4。CasADi也内置了integrator函数来处理更复杂的连续动力学。约束添加opti.subject_to()是添加约束的核心函数。注意我们对每个时间点循环添加约束这可能会产生大量约束条件但CasADi和IPOPT能高效处理稀疏结构。目标函数最小化控制量的平方和即∫ ||a(t)||² dt这通常对应于最小化能量消耗或控制力波动能使轨迹更平滑。初始猜测对于非线性优化问题初始值至关重要。线性插值是一个简单有效的选择。如果问题非凸例如存在障碍物可能需要更复杂的初始化策略甚至使用多起点优化。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, N1)3.4 可视化与结果分析将优化得到的轨迹画出来是验证模型和呈现结果的关键。# 4. 结果可视化 fig plt.figure(figsize(18, 10)) # 1. 3D轨迹图 ax1 fig.add_subplot(231, projection3d) ax1.plot(pos_opt[0, :], pos_opt[1, :], pos_opt[2, :], b-, linewidth2, labelOptimal Trajectory) ax1.scatter(pos_opt[0, 0], pos_opt[1, 0], pos_opt[2, 0], cg, s100, markero, labelStart) ax1.scatter(pos_opt[0, -1], pos_opt[1, -1], pos_opt[2, -1], cr, s100, marker^, labelGoal) for wp in waypoints: _, x, y, z wp ax1.scatter(x, y, z, cm, s80, markers, labelWaypoint 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, :], labelX) ax2.plot(time_axis, pos_opt[1, :], labelY) ax2.plot(time_axis, pos_opt[2, :], labelZ) for wp in waypoints: t_ratio, x, y, z wp t t_ratio * T ax2.axvline(xt, colorgray, linestyle--, alpha0.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, axis0) ax3.plot(time_axis, vel_opt[0, :], labelVx) ax3.plot(time_axis, vel_opt[1, :], labelVy) ax3.plot(time_axis, vel_opt[2, :], labelVz) ax3.plot(time_axis, vel_norm, k--, labelNorm, linewidth2) ax3.axhline(yv_max, colorr, linestyle--, labelfV_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, axis0) # 注意控制量只有N个点时间轴需要调整 time_control np.linspace(0, T, N1)[:-1] dt/2 # 取区间中点代表该区间控制量 ax4.plot(time_control, acc_opt[0, :], labelAx) ax4.plot(time_control, acc_opt[1, :], labelAy) ax4.plot(time_control, acc_opt[2, :], labelAz) ax4.plot(time_control, acc_norm, k--, labelNorm, linewidth2) ax4.axhline(ya_max, colorr, linestyle--, labelfA_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 数学建模论文写作的核心结构论文是将你的思想呈现给评委的唯一载体。结构清晰、逻辑自洽至关重要。一篇完整的数模论文通常包含摘要重中之重用300-500字概括整个工作。必须包含问题重述、你的建模思路、所用方法、主要结果、结论与特色。避免细节突出整体逻辑和创新点。最后写但最先被阅读。问题重述与分析用自己的语言精炼地复述问题并进行分析指出问题的关键点、难点和可能的解决路径。展示你对问题的深刻理解。模型假设与符号说明清晰列出所有假设并说明其合理性。制作一个符号表列出所有主要变量、符号及其含义和单位。模型建立与求解这是论文的躯干。模型准备阐述建模框架如最优控制、数值优化。模型建立详细推导状态方程、约束条件和目标函数。公式要编号解释要清晰。模型求解说明你采用的求解算法如直接转录法IPOPT、离散化方法、以及如何处理特定约束如中间点约束。可以配上算法流程图。参数设置给出模型中所有参数的值及其来源题目给定、合理假设或参考文献。模型求解与结果分析仿真环境说明你的编程语言、工具包和计算平台。基准案例展示一个标准参数下的完整结果就像我们上面代码跑出来的那样配上精心设计的图表和文字分析。分析要指出“由图X可见轨迹平滑穿过所有关键点...”、“表Y显示最大速度未超过约束...”、“这表明模型能有效生成可行轨迹...”。灵敏度分析这是体现模型稳健性和你分析能力的关键改变关键参数如v_max,a_max, 关键点位置总时间T观察目标函数和轨迹形态如何变化。分析其物理意义和工程启示。模型对比/扩展如果时间允许将你的模型与一个简单模型如不考虑能量的最短路径对比或者展示扩展模型如加入障碍物的结果。模型评价与推广客观评价模型的优点如考虑全面、求解高效、结果合理和缺点如假设简化、未考虑XX因素。提出模型的改进方向和在更广领域的应用前景。参考文献规范引用。附录可以放置核心代码不必全部关键片段即可。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)是一个不错的起点。时间管理失控对策三天比赛第一天必须完成问题分析、初步建模和基础代码框架。第二天上午要得到第一个可运行的结果下午进行深入分析和模型改进/扩展。第三天全天用于论文写作、润色和制作图表。编程和建模要并行不要等模型完美了再写代码也不要代码跑通了再开始写论文。论文写作是持续的过程。最后我想强调的是“空中芭蕾”这类题目考察的远不止编程和数学。它考察的是你将模糊需求转化为精确模型的能力、在复杂约束下寻找解决方案的创造力以及用严谨文字和图表展示工作的沟通力。上面的代码和框架提供了一个坚实的起点但真正的亮点来自于你对问题的独特见解和深入分析。多思考“如果...会怎样”并尝试用你的模型去回答你的论文就会脱颖而出。