CasADi 入门实战:符号计算与自动微分加速非线性优化 1. 为什么我要啃 CasADi 这块硬骨头第一次听到 CasADi 这个名字是在一个做机械臂轨迹规划的小圈子里。当时我正被一个六轴机械臂的时间最优轨迹问题折磨得够呛用传统的手写雅可比矩阵加数值优化代码写了八百多行跑一次要等十几秒稍微改一下约束条件就得推倒重来。有个做控制的朋友甩给我一句“你试试 CasADi 吧符号框架加自动微分写起来跟写数学公式差不多。”我当时半信半疑一个开源框架能有多神结果花了一个周末啃完基础文档把原来的问题重写了一遍代码量直接砍到一百多行求解速度还快了一个数量级。从那以后CasADi 就成了我工具箱里常驻的核心组件。CasADi 到底是什么简单说它是一个用于非线性优化和最优控制的开源框架核心能力是符号计算加自动微分。你可以把它理解成一个“数学公式翻译器加求导器”你用接近数学表达式的语法把问题写出来它自动帮你算出雅可比矩阵、海森矩阵然后对接各种数值求解器比如 IPOPT、SNOPT、qpoases去求解。它解决的问题非常明确——当你需要处理带约束的非线性优化问题尤其是涉及微分方程、动态系统、机器人运动学动力学这类场景时手写导数既容易出错又极其耗时CasADi 就是来把你从这种重复劳动里解放出来的。这篇文章适合谁看如果你正在做机器人控制、轨迹优化、模型预测控制MPC、系统辨识、或者任何需要反复求解非线性规划NLP的活儿而且你已经被手动求导或者数值差分折磨过那这篇内容就是为你准备的。如果你只是听说过 CasADi 但一直没敢下手或者装完环境跑了个 demo 就不知道下一步该干嘛那正好我把从零上手到跑通第一个实际问题的完整路径都梳理一遍。我不会只给你贴官方文档的翻译而是把我在实际项目里踩过的坑、绕过的弯、总结出来的高效用法都摊开来讲。2. CasADi 的核心设计思路拆解2.1 符号框架与数值求解的边界在哪里很多人第一次接触 CasADi 会有一个误解以为它是个“求解器”。其实不是。CasADi 本身不做数值求解它是一个建模层和微分层。真正的求解工作是由它背后对接的求解器完成的比如 IPOPT 负责非线性规划qpoases 负责二次规划SUNDIALS 负责微分方程积分。CasADi 的职责是让你用一种统一的符号语言描述问题自动生成求解器需要的所有导数信息然后把问题打包成求解器能吃的格式递过去。这个设计思路的好处非常明显。传统做法里你写一个优化问题目标函数和约束的梯度、雅可比、海森矩阵都得自己推公式、自己写代码。一旦问题结构变了比如加了一个新约束所有导数都得重新推。CasADi 用符号表达式记录你的计算图然后通过自动微分具体来说是反向模式 AD 和前向模式 AD 的结合在需要的时候精确计算任意阶导数。注意是精确计算不是数值差分近似。数值差分不仅慢而且步长选不好精度直接崩掉。自动微分没有这个问题它利用链式法则在计算图上一层一层传播导数精度是机器精度级别。那什么时候该用 CasADi什么时候不该用我的经验是如果你的问题可以用矩阵运算和基本数学函数sin、cos、exp、log、sqrt 等表达出来而且你需要反复求导或者反复求解那 CasADi 就是利器。反过来如果你的问题本身很简单比如一个无约束的二次函数最小化手写解析解或者用现成的库更省事。CasADi 的价值在复杂问题和反复迭代的场景里才真正体现出来。2.2 自动微分为什么比手写导数靠谱我举个实际例子你就明白了。假设你要做一个两连杆机械臂的逆动力学计算末端执行器的位置是关节角度的函数p l1*cos(q1) l2*cos(q1q2)你要算末端位置对 q1 和 q2 的偏导数手写的话就是dp/dq1 -l1*sin(q1) - l2*sin(q1q2) dp/dq2 -l2*sin(q1q2)这还只是两个关节。如果是六轴机械臂末端位置对六个关节角的偏导数手写雅可比矩阵就是一场噩梦而且极容易在某个符号上出错。更麻烦的是如果你后面还要算二阶导数比如做基于梯度的优化需要海森矩阵手写几乎不可能保证正确。CasADi 的做法是你只管写p l1*cos(q1) l2*cos(q1q2)然后调用jacobian(p, [q1, q2])它自动给你返回符号形式的雅可比表达式。你还可以进一步对这个雅可比再求导得到二阶导数。整个过程不需要你手动推任何公式。我在实际项目里做过对比一个包含 12 个状态变量、6 个控制输入的 MPC 问题手写导数的版本调试了整整两天才跑通而且每次改模型都要重新检查导数用 CasADi 重写后模型改完直接跑导数自动更新再也没有因为导数错误导致求解器不收敛的情况。注意自动微分虽然精确但符号表达式如果写得过于复杂比如大量嵌套的条件判断会导致计算图膨胀求导效率下降。我的经验是尽量用向量化操作避免在符号表达式里写循环。2.3 与同类工具的横向对比市面上做优化建模的工具不止 CasADi 一个。我列几个常见的对比一下方便你判断什么场景该选什么。工具核心特点适用场景学习曲线CasADi符号建模 自动微分 多求解器接口最优控制、MPC、机器人中等YALMIPMATLAB 建模层侧重凸优化凸优化、LMI较低MATLAB 用户CVXPYPython 凸优化建模凸优化、投资组合较低Gurobi/CPLEX商业求解器自带建模接口线性/整数规划中等PyomoPython 通用优化建模通用数学规划中等CasADi 的差异化优势在于非凸问题和微分方程约束。YALMIP 和 CVXPY 主要面向凸优化虽然也能处理一些非凸问题但能力和 CasADi 不在一个量级。Gurobi 和 CPLEX 是求解器建模能力相对有限。Pyomo 更偏向于传统数学规划对动态系统的支持不如 CasADi 自然。所以如果你的问题涉及微分方程、动态系统、或者需要反复求高阶导数CasADi 基本是首选。3. 环境搭建与第一个可运行示例3.1 安装路径选择与依赖处理CasADi 的安装方式主要有三种pip 安装、conda 安装、源码编译。我强烈建议新手走pip 安装一条命令搞定pip install casadi如果你用 conda也可以conda install -c conda-forge casadi源码编译只推荐给需要自定义求解器接口或者做二次开发的人普通用户没必要折腾。编译过程涉及 CMake、IPOPT 源码、HSL 库等一系列依赖我第一次编译花了整整一个下午中间因为 HSL 库的许可证问题卡了很久。除非你有明确需求否则 pip 版本足够用。安装完成后验证一下import casadi as ca print(ca.__version__)如果输出版本号说明安装成功。这里有个小坑CasADi 的 Python 包名是casadi但导入时通常简写为ca这是社区惯例官方文档也这么用。提示如果你在 Windows 上安装建议用 Python 3.8 到 3.11 之间的版本太新的 Python 版本可能还没有预编译的 wheel 包。我在 Python 3.12 刚发布时试过pip 找不到对应 wheel只能降级。3.2 符号变量定义与基本运算CasADi 的核心对象是MX和SX。这两个都是符号类型区别在于SX标量符号适合小规模问题展开后计算效率高。MX矩阵符号适合大规模问题支持稀疏性但单个元素的操作开销比 SX 大。新手先用 SX 就行简单直观。定义一个符号变量import casadi as ca x ca.SX.sym(x) y ca.SX.sym(y)这就定义了两个标量符号。你也可以定义向量v ca.SX.sym(v, 3) # 3维向量 M ca.SX.sym(M, 2, 2) # 2x2矩阵定义好符号后就可以像写普通数学表达式一样构造函数f x**2 ca.sin(y) x*y这个f现在是一个符号表达式不是数值。你可以对它求导df_dx ca.jacobian(f, x) df_dy ca.jacobian(f, y)ca.jacobian返回的也是符号表达式。如果你想把它变成可调用的函数用ca.Functionf_func ca.Function(f, [x, y], [f]) print(f_func(2, 3)) # 输出数值结果这个ca.Function是 CasADi 里非常核心的概念后面做优化、积分、求解都要用到它。你可以把它理解成一个“编译后的数学函数”输入数值输出数值内部的计算图已经优化过了。3.3 一个完整的非线性规划求解实例光说不练假把式。我们来看一个完整的例子求解一个带约束的非线性优化问题。问题描述最小化(x1 - 3)^2 (x2 - 2)^2约束条件是x1^2 x2^2 5x1 0x2 0。用 CasADi 写出来是这样的import casadi as ca # 定义符号变量 x1 ca.SX.sym(x1) x2 ca.SX.sym(x2) x ca.vertcat(x1, x2) # 拼成向量 # 目标函数 obj (x1 - 3)**2 (x2 - 2)**2 # 约束条件 g x1**2 x2**2 # 不等式约束 g 5 lbg -ca.inf ubg 5 # 变量上下界 lbx [0, 0] ubx [ca.inf, ca.inf] # 构造 NLP 问题 nlp {x: x, f: obj, g: g} solver ca.nlpsol(solver, ipopt, nlp) # 初始猜测 x0 [1, 1] # 求解 sol solver(x0x0, lbxlbx, ubxubx, lbglbg, ubgubg) print(sol[x])跑一下输出应该是接近[2.236, 0]或者[1.8, 1.6]附近的解具体取决于初始点。这个例子里我们只写了目标函数和约束的表达式没有写任何导数。IPOPT 求解器需要的梯度、雅可比、海森矩阵CasADi 全部自动生成。这就是它最核心的价值。注意ca.nlpsol的第一个参数是求解器名字第二个参数是求解器类型。IPOPT 是开源里最常用的非线性求解器适合中小规模问题。如果你的问题有大量整数变量IPOPT 就不合适了得换别的工具。4. 实操过程中容易踩的坑与排查技巧4.1 求解器不收敛的常见原因第一次用 CasADi 跑优化十有八九会遇到求解器不收敛的情况。IPOPT 会返回一个状态码比如success、max_iter、infeasible之类的。我整理了一个排查表按优先级从高到低检查现象可能原因排查方法迭代次数达到上限初始点太差换几个不同的初始点试试约束不可行约束条件互相矛盾检查约束边界是否合理梯度计算异常表达式里有除零或 log 负数检查符号表达式定义域求解速度极慢问题规模太大或表达式太复杂用 MX 替代 SX启用稀疏性返回 NaN数值溢出检查是否有指数爆炸项我印象最深的一次是做一个轨迹优化求解器一直返回infeasible。查了半天发现是约束里有一个sqrt函数而某个中间变量的值在迭代过程中变成了负数导致sqrt返回 NaN整个求解崩溃。解决办法是在sqrt外面套一个ca.fmax(0, x)保证输入非负。这种问题在纯数值代码里也会遇到但在符号框架里更容易被忽略因为表达式是“一次性写好”的运行时才暴露问题。4.2 符号表达式膨胀与性能优化CasADi 的符号表达式如果写得不好会变得极其庞大。比如你在符号层面写了一个循环每次迭代都往表达式里追加一项最后这个表达式可能有几千个节点。求导的时候计算图会爆炸内存和时间都吃不消。我的经验是尽量用向量化操作避免在符号层面写循环。比如你要计算一个求和# 不好的写法 total 0 for i in range(100): total x[i]**2 # 好的写法 total ca.dot(x, x)ca.dot是 CasADi 内置的向量点积生成的符号表达式简洁得多。类似的内置函数还有ca.sum1、ca.sum2、ca.mtimes等。另外如果问题规模确实很大考虑用MX代替SX。MX支持稀疏矩阵运算不会把每个元素都展开成标量节点内存占用和求导效率都会好很多。还有一个技巧是把重复使用的子表达式提取出来。比如# 重复计算 f1 ca.sin(x) ca.cos(x) f2 ca.sin(x) * ca.cos(x) # 提取公共项 s ca.sin(x) c ca.cos(x) f1 s c f2 s * c这样 CasADi 在构建计算图时会共享子节点减少冗余。4.3 与 NumPy 的互操作细节CasADi 的符号类型和 NumPy 数组之间的转换是新手经常卡住的地方。几个关键点ca.DM是 CasADi 的数值矩阵类型可以和 NumPy 数组互相转换。ca.DM(np_array)把 NumPy 数组转成 CasADi 矩阵。np.array(ca_dm)把 CasADi 矩阵转回 NumPy 数组。符号表达式SX/MX不能直接转 NumPy必须先变成ca.Function再调用。我踩过的一个坑是在定义优化问题时把 NumPy 数组直接传给了ca.SX.sym的维度参数结果报错。正确的做法是用整数指定维度或者用ca.DM包装数值数据。另外CasADi 的矩阵索引是从 0 开始的和 NumPy 一致但切片语法略有不同建议查一下官方文档的索引规则。提示如果你在 Jupyter Notebook 里用 CasADi记得把符号表达式用ca.Function包装后再打印否则输出的是巨大的符号树Notebook 会卡死。我吃过这个亏一个不小心打印了一个未包装的表达式浏览器直接无响应。5. 从学习到实战的进阶路径5.1 掌握 CasADi 需要补哪些数学基础CasADi 本身不难学但它背后的数学概念需要你有一定基础。我列几个关键点非线性规划理解目标函数、约束、拉格朗日乘子、KKT 条件。不需要你会手推但至少知道求解器在干什么。自动微分理解链式法则、前向模式、反向模式。知道为什么自动微分比数值差分准。常微分方程如果你要做最优控制得知道 ODE 的基本形式和数值积分方法如 RK4。稀疏矩阵大规模问题里稀疏性决定性能知道雅可比矩阵的稀疏结构怎么影响求解效率。这些不需要你一次性全学会边用边补就行。我的路径是先跑通官方示例然后拿一个自己手头的小问题练手遇到不懂的概念再回头查资料。这样学得最快也最扎实。5.2 官方示例与文档的高效使用方法CasADi 的官方文档质量很高但内容很多新手容易迷失。我的建议是先看User Guide的前三章把 SX、MX、Function 这三个概念搞清楚。然后直接跳到Examples页面找和你问题最接近的示例复制下来跑通。跑通后逐行修改参数观察输出变化理解每个参数的作用。最后再回头补理论细节。官方示例里有一个optimal_control的目录里面包含了 MPC、轨迹优化、微分方程约束优化的完整代码。这些示例是我见过的最好的学习材料比任何教程都管用。我当初就是把race_car那个示例反复改了十几遍才真正理解 CasADi 在最优控制里的用法。5.3 一个真实项目的简化复盘最后分享一个我用 CasADi 做的实际项目片段。项目背景是一个三自由度的平面机械臂需要做点到点的最优轨迹规划约束包括关节角度限制、关节速度限制、以及末端执行器不能碰到一个圆形障碍物。核心代码结构是这样的import casadi as ca # 状态变量q1, q2, q3, dq1, dq2, dq3 q ca.SX.sym(q, 3) dq ca.SX.sym(dq, 3) x ca.vertcat(q, dq) # 控制变量关节力矩 u ca.SX.sym(u, 3) # 动力学方程简化版 M ca.SX.zeros(3, 3) # 质量矩阵实际项目里是 q 的函数 # ... 填充 M ... ddq ca.solve(M, u) # 加速度 # 状态方程 xdot ca.vertcat(dq, ddq) f ca.Function(f, [x, u], [xdot]) # 障碍物约束 p_end forward_kinematics(q) # 末端位置 obstacle_center ca.DM([0.5, 0.5]) dist ca.norm_2(p_end - obstacle_center) g_obs dist # 约束dist 0.2 # 构造最优控制问题 # ... 离散化、构造 NLP、调用 IPOPT ...这个项目里CasADi 帮我省掉了所有动力学方程的符号推导和导数计算。质量矩阵 M 是 q 的复杂函数手写导数几乎不可能但 CasADi 自动处理了。整个开发周期从预计的两周压缩到三天而且代码可读性极高后续改约束条件只需要改几行。注意动力学方程里的ca.solve(M, u)在 M 接近奇异时会出问题。实际项目里我加了一个正则化项M 1e-6*ca.SX.eye(3)避免数值不稳定。这种细节官方文档不会写但实际用的时候必须考虑。6. 一些零散但重要的经验补充6.1 版本兼容性与更新策略CasADi 的版本更新比较频繁不同版本之间的 API 偶尔会有变化。我的建议是生产项目锁定版本学习阶段用最新版。锁定版本可以避免因为 API 变动导致代码突然跑不了。我一般会在requirements.txt里写死版本号比如casadi3.6.3。学习阶段用最新版可以体验到新功能和性能改进。另外CasADi 的 Python 接口和 C 接口在功能上基本对齐但有些高级特性比如自定义求解器接口只在 C 里提供。如果你只是做建模和求解Python 版完全够用。6.2 调试符号表达式的实用技巧调试符号表达式比调试数值代码麻烦因为你看不到中间值。我的几个常用技巧用ca.Function把中间表达式包装成函数传入具体数值打印输出验证逻辑。用ca.gradient和ca.jacobian检查导数表达式是否符合预期。对于复杂表达式分段验证先验证一小部分再逐步拼接。用ca.print_expression或者直接print符号表达式但注意表达式太大时会刷屏。我遇到过一个 bug符号表达式里用了ca.if_else但两个分支的维度不一致导致运行时维度错误。这种问题在数值代码里很容易发现但在符号层面需要仔细检查每个分支的输出形状。6.3 社区资源与持续学习建议CasADi 的社区不算特别大但质量很高。官方论坛里有很多深度讨论尤其是关于求解器调参和性能优化的帖子。另外GitHub 上的 issue 区也值得翻很多常见问题别人已经踩过坑了。如果你想深入我建议读一读 CasADi 作者的那篇论文Andersson et al., 2019里面详细讲了框架的设计哲学和自动微分的实现细节。读完之后你对 CasADi 的理解会从“会用”上升到“知道为什么这么设计”。最后再分享一个小技巧CasADi 的ca.Function支持jit选项可以启用即时编译加速。对于反复调用的函数开启jitTrue能显著提升性能。不过 JIT 编译第一次调用会慢一些适合长期运行的任务。这个内容后续还可以这样扩展把 CasADi 和深度学习框架结合用自动微分做端到端的优化控制或者用 CasADi 做模型预测控制对接实际硬件做实时控制。这些方向我都在探索后面有机会再单独写。