SymPy Physics/Mechanics 线性化实战指南用 Linearizer 在平衡点附近求解 Kane 与 Lagrange 方程的线性化运动方程【免费下载链接】sympyA computer algebra system written in pure Python项目地址: https://gitcode.com/GitHub_Trending/sy/sympysympy.physics.mechanics提供了围绕工作点operating point也称 trim condition对生成的运动方程EOM, equations of motion进行线性化的完整工具链。本文基于 SymPy 官方文档 linearize.rst 展开系统讲解线性化的数学背景、M/A/B与A/B两种输出形式、KanesMethod与LagrangesMethod两条使用路径并结合 linearize.py 源码剖析置换矩阵、系数矩阵与线性求解器的内部机制最后给出常见的性能与nan/zoo/oo排障技巧。读完本文你可以直接复制文中的单摆示例对任意受约束多体系统进行稳定性分析或控制器设计所需的状态空间线性化。线性化的数学背景广义系统形式在sympy.physics.mechanics中假设所有系统都可以表示为如下广义形式五个动力/运动方程组加三个约束方程组$$ \begin{aligned} f_{c}(q, t) 0_{l \times 1}\ f_{v}(q, u, t) 0_{m \times 1}\ f_{a}(q, \dot{q}, u, \dot{u}, t) 0_{m \times 1}\ f_{0}(q, \dot{q}, t) f_{1}(q, u, t) 0_{n \times 1}\ f_{2}(q, u, \dot{u}, t) f_{3}(q, \dot{q}, u, r, t) f_{4}(q, \lambda, t) 0_{(o-mk) \times 1} \end{aligned} $$其中各向量的维度为$$ \begin{aligned} q, \dot{q} \in \mathbb{R}^n\ u, \dot{u} \in \mathbb{R}^o\ r \in \mathbb{R}^s\ \lambda \in \mathbb{R}^k \end{aligned} $$各符号的含义符号含义$f_c$构型约束方程configuration constraints$f_v$速度约束方程velocity constraints$f_a$加速度约束方程acceleration constraints$f_0$、$f_1$运动学微分方程$f_2$、$f_3$、$f_4$动力学微分方程$q$、$\dot{q}$广义坐标及其导数$u$、$\dot{u}$广义速度generalized speeds及其导数$r$系统输入$\lambda$Lagrange 乘子这个广义形式由 linearize.py 中的Linearizer类持有并由它完成实际的线性化工作。KanesMethod和LagrangesMethod对象都提供to_linearizer类方法把各自的方法框架Kane 方程或 Lagrange 方程强制转换coerce为上述广义形式。关于从属坐标与从属速度的说明如果待线性化的系统包含约束方程则并非所有广义坐标都是独立的例如 $q_1$ 可能依赖于 $q_2$。设有 $l$ 个构型约束和 $m$ 个速度约束则存在 $l$ 个从属坐标dependent coordinates和 $m$ 个从属速度dependent speeds。从源码结构看这一划分在Linearizer.__init__中对应q_i/q_d独立/从属坐标与u_i/u_d独立/从属速度四个向量维度元组(l, m, n, o, s, k)会在构造时一次性推导出来见 linearize.py 的self._dims。一般地你可以任选哪些坐标/速度作为从属变量但实践中某些选择会在特定工作点引发数值奇异。如何系统地决定从属变量划分超出本文范围文档建议参考 Blajer (1994) 的坐标划分方法。线性化方程的两种输出形式sympy.physics.mechanics提供两种线性化 EOM 的输出形式。形式一M、A、B隐式形式默认在这种形式中力矩阵被线性化为两个独立矩阵 $A$ 和 $B$这是默认的线性化输出。方程为$$ M \begin{bmatrix} \delta \dot{q} \ \delta \dot{u} \ \delta \lambda \end{bmatrix} A \begin{bmatrix} \delta q_i \ \delta u_i \end{bmatrix} B \begin{bmatrix} \delta r \end{bmatrix} $$其中$$ \begin{aligned} M \in \mathbb{R}^{(nok) \times (nok)}\ A \in \mathbb{R}^{(nok) \times (n-lo-m)}\ B \in \mathbb{R}^{(nok) \times s} \end{aligned} $$注意 $q_i$ 和 $u_i$ 仅指独立坐标与速度而 $q$ 和 $u$ 同时包含独立与从属变量。这种形式保留了解线性方程组的自由度——你可以先符号化得到 M、A、B之后再代入数值工作点手工求解是大型符号系统的推荐做法。形式二A、B显式一阶状态空间形式这种形式把线性化 EOM 化为仅含独立坐标与速度的显式一阶形式常用于稳定性分析或控制理论$$ \begin{bmatrix} \delta \dot{q_i} \ \delta \delta \dot{u_i} \end{bmatrix} A \begin{bmatrix} \delta q_i \ \delta u_i \end{bmatrix} B \begin{bmatrix} \delta r \end{bmatrix} $$其中$$ \begin{aligned} A \in \mathbb{R}^{(n-lo-m) \times (n-lo-m)}\ B \in \mathbb{R}^{(n-lo-m) \times s} \end{aligned} $$使用该形式时在linearize方法中设置A_and_BTrue即可。对 Kane 方程线性化单摆示例下面用一个简单单摆系统演示完整流程。首先用KanesMethod建立系统并生成F_r与F_r^* from sympy import symbols, Matrix from sympy.physics.mechanics import * q1 dynamicsymbols(q1) # 摆的角度 u1 dynamicsymbols(u1) # 角速度 q1d dynamicsymbols(q1, 1) L, m, t, g symbols(L, m, t, g) # 构造世界坐标系 N ReferenceFrame(N) pN Point(N*) pN.set_vel(N, 0) # A.x 沿摆杆方向 A N.orientnew(A, axis, [q1, N.z]) A.set_ang_vel(N, u1*N.z) # 相对原点 N* 定位点 P P pN.locatenew(P, L*A.x) vel_P P.v2pt_theory(pN, N, A) pP Particle(pP, P, m) # 运动学微分方程 kde Matrix([q1d - u1]) # 作用在 P 点的力此处文档示例为水平力 R m*g*N.x # 用 Kane 法求解运动方程 KM KanesMethod(N, q_ind[q1], u_ind[u1], kd_eqskde) fr, frstar KM.kanes_equations([pP], [(P, R)])方式一直接使用 Linearizer 类通过to_linearizer类方法创建 linearizer 对象。该方法会把KanesMethod对象中的表示转换为上文描述的广义形式。由于独立/从属坐标与速度在创建KanesMethod时已经指定这里无需再次指定 linearizer KM.to_linearizer()然后调用Linearizer对象的linearize方法生成线性化 EOM M, A, B linearizer.linearize() M Matrix([ [1, 0], [0, -L**2*m]]) A Matrix([ [ 0, 1], [L*g*m*cos(q1(t)), 0]]) B Matrix(0, 0, [])也可以通过指定A_and_BTrue生成 A、B 形式 A, B linearizer.linearize(A_and_BTrue) A Matrix([ [ 0, 1], [-g*cos(q1(t))/L, 0]]) B Matrix(0, 0, [])工作点可以指定为字典或字典的可迭代对象。系统会在返回矩阵之前先在该指定点上求值 op_point {q1: 0, u1: 0} A_op, B_op linearizer.linearize(A_and_BTrue, op_pointop_point) A_op Matrix([ [ 0, 1], [-g/L, 0]])同样的效果也可以通过把msubs应用到未指定op_point时生成的矩阵上来实现 assert msubs(A, op_point) A_op有时返回的矩阵未必处于最简形式。你可以事后再做化简或者把Linearizer对象的simplify关键字参数设为True让它在内部执行化简。方式二使用 linearize 类方法KanesMethod类还提供了linearize方法作为便捷封装它内部调用to_linearizer执行线性化并返回结果。上文linearize方法中可用的所有关键字参数在这里同样可用 A, B, inp_vec KM.linearize(A_and_BTrue, op_pointop_point) A Matrix([ [ 0, 1], [-g/L, 0]])从源码看kane.py 中的KanesMethod.linearize实现为三行核心逻辑先self.to_linearizer(linear_solverlinear_solver)再linearizer.linearize(**kwargs)最后把输入向量linearizer.r追加到结果元组中见 kane.py。额外的输出inp_vec是一个包含所有未被广义坐标或速度向量收录的dynamicsymbols的向量。这些符号被视为系统的输入构成背景部分描述的 $r$ 向量。本例中没有输入所以向量为空 inp_vec Matrix(0, 0, [])从源码结构看r向量的构造发生在to_linearizer内部KanesMethod用find_dynamicsymbols在动力学驱动项f_3中搜索不属于q、qdot、u、udot、辅助速度及其导数的所有动态符号并用default_sort_key按正则序排序从而得到确定性的输入向量见 kane.py。注意早期 SymPy 版本中的KanesMethod线性化方法已被弃用将在未来版本移除。当前仓库中KanesMethod.linearize保留了new_method关键字参数用于向后兼容从源码看它已标注为 Deprecated, does nothing and will be removed见 kane.py新代码无需再显式传new_methodTrue。对 Lagrange 方程线性化Lagrange 方程的线性化流程与 Kane 方程基本相同。同样以单摆为例先把 A 坐标系与 P 点用q1d而非u1重新定义 # 用 q1d 而非 u1 重新定义 A 和 P A N.orientnew(A, axis, [q1, N.z]) A.set_ang_vel(N, q1d*N.z) P pN.locatenew(P, L*A.x) vel_P P.v2pt_theory(pN, N, A) pP Particle(pP, P, m) # 用 Lagrange 法求解运动方程 Lag Lagrangian(N, pP) LM LagrangesMethod(Lag, [q1], forcelist[(P, R)], frameN) lag_eqs LM.form_lagranges_equations()方式一直接使用 Linearizer 类LagrangesMethod对象同样可以用to_linearizer类方法生成Linearizer对象。与KanesMethod的唯一区别在于LagrangesMethod对象内部并没有预先指定独立/从属坐标与速度这些必须在调用to_linearizer时显式给出。本例没有从属坐标和速度如果有它们会包含在q_dep和qd_dep关键字参数中 linearizer LM.to_linearizer(q_ind[q1], qd_ind[q1d])此后流程与KanesMethod示例完全一致 A, B linearizer.linearize(A_and_BTrue, op_pointop_point) A Matrix([ [ 0, 1], [-g/L, 0]])方式二使用 linearize 类方法与KanesMethod类似LagrangesMethod类也提供linearize方法作为封装内部调用to_linearizer、执行线性化并返回结果。同样地独立/从属坐标与速度必须在调用时给出 A, B, inp_vec LM.linearize(q_ind[q1], qd_ind[q1d], A_and_BTrue, op_pointop_point) A Matrix([ [ 0, 1], [-g/L, 0]])从源码结构看lagrange.py 中的LagrangesMethod.to_linearizer会把 Lagrange 方程的各部分映射到广义形式f_0 u、f_1 -u运动学方程、f_2 _term1质量矩阵项、f_3 -(_term2 _term4)驱动力项、f_4 -_term3Lagrange 乘子项对应 $f_4$并把hol_coneqs与coneqs分别映射为 $f_c$ 与 $f_v$、对其求时间导数得到 $f_a$见 lagrange.py。该函数还会严格校验从属变量划分q_dep数量必须等于构型约束数、qd_dep数量必须等于速度约束数且q_ind q_dep必须恰好是q的一个划分否则抛出ValueError见 lagrange.py。另外带 Lagrange 乘子的 Lagrange 系统还可用solve_multipliers在指定工作点求解乘子值这在 test_linearize.py 的 nonminimal 单摆测试中被用于补全工作点字典属于约束系统线性化工作流的一部分。Linearizer 的内部机制置换矩阵、系数矩阵与线性求解器上面的 API 背后Linearizer做了几件重要的事理解它们有助于在大型系统上做出正确决策。延迟初始化_setup构造Linearizer时只保存广义方程向量与各变量向量置换矩阵、块矩阵、系数矩阵的计算被推迟到首次调用linearize时才执行_setup_done标志位以提高构造速度见 linearize.py。置换矩阵_form_permutation_matrices利用permutation_matrix构造Pq、Pu以及组合置换矩阵perm_mat满足 $[q_{ind}, u_{ind}]^T P\cdot[q, u]^T$。从属坐标/速度的排序约定是独立变量在前、从属变量在末尾P_qi/P_qd等子块按l、m切分见 linearize.py。系数矩阵 C_0、C_1、C_2_form_coefficient_matrices处理约束的影响。有构型约束$l 0$时 $C_0 (I - P_{qd}\cdot\mathrm{solver}(f_{c,jac},P_{qd}, f_{c,jac}))P_{qi}$无约束时 $C_0 I$。有运动约束$m 0$时类似地构造 $C_1$、$C_2$否则 $C_1 0$、$C_2 I$见 linearize.py。这一步正是正确处理约束的落点——约束导致从属变量与独立变量耦合一阶 Taylor 展开不能简单对右端求 Jacobian。块矩阵拼装_form_block_matrices对 $f_0, f_1, f_2, f_3, f_4$ 及各约束方程求 Jacobian按存在性条件拼出 $M$、$A$、$B$ 的各个子块如M_qq f_0.jacobian(qd)、A_qu -f_1.jacobian(u)等见 linearize.py。linearize方法再按注释中的块结构把这些子块row_join成完整矩阵。可插拔线性求解器linear_solver参数既可以是MatrixBase.solve接受的方法名字符串默认LU对应A.LUsolve(b)也可以是形如x f(A, b)的可调用对象。_parse_linear_solver负责解析见 functions.py。文档特别提示LUsolve计算快但常常遇到零除而产生nan结果——这正是排障章节第二个问题的底层原因之一。linearize方法的关键字参数汇总参数默认值说明op_pointNone工作点字典或字典的可迭代对象值可为数值或符号符号被替换得越多运行越快A_and_BFalseFalse返回(M, A, B)True返回(A, B)simplifyFalse是否在返回前对矩阵做符号化简对大表达式可能非常耗时潜在问题与排障虽然Linearizer类理论上能线性化所有系统实践中仍会遇到一些问题。问题一符号化线性化且 A_and_BTrue 时很慢最可能的原因是符号化求解大型线性方程组本身是昂贵操作。指定工作点可以缩小表达式规模、加速求解。如果确实需要纯符号解例如之后要批量代入许多工作点可以绕过A_and_BTrue的内部符号求解先用A_and_BFalse取回 M、A、B代入工作点后再手工求解 M, A, B linearizer.linearize() M_op msubs(M, op_point) A_op msubs(A, op_point) perm_mat linearizer.perm_mat A_lin perm_mat.T * M_op.LUsolve(A_op) A_lin Matrix([ [ 0, 1], [-g/L, 0]])这与源码注释给出的公式一致$A P^T M^{-1} A$、$B P^T M^{-1} B$其中 $P \texttt{Linearizer.perm_mat}$见 linearize.py。求解前A、M中符号越少求解越快对大型表达式把转换到 A、B 形式推迟到大部分符号已代入数值之后收益明显。问题二线性化结果出现 nan、zoo 或 oo有两种可能的成因。成因一优先检查从属坐标选择在特定工作点导致奇异。某些坐标划分在特定工作点会产生奇异。系统性地避开它需要坐标划分理论超出本文范围参见 Blajer 1994。成因二代入工作点前矩阵未充分化简。一个典型的演示 from sympy import sin, tan expr sin(q1)/tan(q1) op_point {q1: 0} expr.subs(op_point) nan如果先化简再代入就能得到正确值 expr.simplify().subs(op_point) 1目前尚未找到完全自动避免该问题的通用办法。对规模合理的表达式用msubs并设smartTrue会启用一种规避此类条件的算法但对大型表达式它极其耗时 msubs(expr, op_point, smartTrue) 1从源码结构看msubs是physics.mechanics专用的自研替换函数见 functions.py它只遍历表达式树一次跳过Derivative内部项并且支持传入多个替换字典smartTrue时会额外检查直接代入会产生nan、但化简后有效的条件。验证与进阶测试套件中的非最小坐标示例上述单摆示例虽简单但没有包含从属坐标或从属速度。仓库测试 test_linearize.py 提供了可直接运行的更复杂示例nonminimal 单摆Kane 法用直角坐标q1, q2摆锤的 N.x、N.y 坐标与速度u1, u2建模显式给出构型约束f_c |P.pos_from(pN)| - L与速度约束f_v P.vel(N).express(A).dot(A.x)通过q_dependent[q1], u_dependent[u1]指定从属变量在直立工作点线性化后得到A [[0, 1], [-9.8/L, 0]]见 test_linearize.py。该测试同时验证了linear_solverGJ与自定义 callable 求解器三种方式的等价性。nonminimal 单摆Lagrange 法用hol_coneqsf_c建立约束先用solve_multipliers求出 Lagrange 乘子的工作点值并入op_point再线性化验证见 test_linearize.py。滚动圆盘Kane 法含加速度约束6 个广义坐标、6 个广义速度的非完整约束系统在直立滚动名义点线性化后验证 A 矩阵并检查临界速度处特征值全为零见 test_linearize.py。这个用例展示了op_point接收多字典列表[q_op, u_op, qd_op, ud_op]的用法。此外文档还指引读者参考 Tutorials 页中 Nonminimal Coordinates Pendulum 一节其中用 Kane 与 Lagrange 两种方法对同一单摆分别带从属坐标做了线性化可作为本指南示例之外的延伸阅读。小结SymPy 的线性化路径可以概括为一条清晰的调用链KanesMethod/LagrangesMethod生成 EOM →to_linearizer强制转换为 $f_c/f_v/f_a/f_0\sim f_4$ 广义形式自动完成约束感知→Linearizer.linearize做一阶 Taylor 展开并按块矩阵拼装。无约束时它退化为右端对 $q$、$u$ 的 Jacobian有约束时通过置换矩阵与 $C_0, C_1, C_2$ 系数矩阵正确处理从属变量。实践要点大型符号系统优先用M/A/B形式加op_point数值代入再手工求解出现nan时先怀疑坐标划分其次用msubs(..., smartTrue)或化简后再代入。相关实现与测试分别位于 linearize.py、kane.py、lagrange.py 与 test_linearize.py可按需深入。【免费下载链接】sympyA computer algebra system written in pure Python项目地址: https://gitcode.com/GitHub_Trending/sy/sympy创作声明：本文部分内容由AI辅助生成（AIGC），仅供参考 SymPy Physics/Mechanics 线性化实战指南：用 Linearizer 在平衡点附近求解 Kane 与 Lagrange 方程的线性化运动方程
SymPy Physics/Mechanics 线性化实战指南用 Linearizer 在平衡点附近求解 Kane 与 Lagrange 方程的线性化运动方程【免费下载链接】sympyA computer algebra system written in pure Python项目地址: https://gitcode.com/GitHub_Trending/sy/sympysympy.physics.mechanics提供了围绕工作点operating point也称 trim condition对生成的运动方程EOM, equations of motion进行线性化的完整工具链。本文基于 SymPy 官方文档 linearize.rst 展开系统讲解线性化的数学背景、M/A/B与A/B两种输出形式、KanesMethod与LagrangesMethod两条使用路径并结合 linearize.py 源码剖析置换矩阵、系数矩阵与线性求解器的内部机制最后给出常见的性能与nan/zoo/oo排障技巧。读完本文你可以直接复制文中的单摆示例对任意受约束多体系统进行稳定性分析或控制器设计所需的状态空间线性化。线性化的数学背景广义系统形式在sympy.physics.mechanics中假设所有系统都可以表示为如下广义形式五个动力/运动方程组加三个约束方程组$$ \begin{aligned} f_{c}(q, t) 0_{l \times 1}\ f_{v}(q, u, t) 0_{m \times 1}\ f_{a}(q, \dot{q}, u, \dot{u}, t) 0_{m \times 1}\ f_{0}(q, \dot{q}, t) f_{1}(q, u, t) 0_{n \times 1}\ f_{2}(q, u, \dot{u}, t) f_{3}(q, \dot{q}, u, r, t) f_{4}(q, \lambda, t) 0_{(o-mk) \times 1} \end{aligned} $$其中各向量的维度为$$ \begin{aligned} q, \dot{q} \in \mathbb{R}^n\ u, \dot{u} \in \mathbb{R}^o\ r \in \mathbb{R}^s\ \lambda \in \mathbb{R}^k \end{aligned} $$各符号的含义符号含义$f_c$构型约束方程configuration constraints$f_v$速度约束方程velocity constraints$f_a$加速度约束方程acceleration constraints$f_0$、$f_1$运动学微分方程$f_2$、$f_3$、$f_4$动力学微分方程$q$、$\dot{q}$广义坐标及其导数$u$、$\dot{u}$广义速度generalized speeds及其导数$r$系统输入$\lambda$Lagrange 乘子这个广义形式由 linearize.py 中的Linearizer类持有并由它完成实际的线性化工作。KanesMethod和LagrangesMethod对象都提供to_linearizer类方法把各自的方法框架Kane 方程或 Lagrange 方程强制转换coerce为上述广义形式。关于从属坐标与从属速度的说明如果待线性化的系统包含约束方程则并非所有广义坐标都是独立的例如 $q_1$ 可能依赖于 $q_2$。设有 $l$ 个构型约束和 $m$ 个速度约束则存在 $l$ 个从属坐标dependent coordinates和 $m$ 个从属速度dependent speeds。从源码结构看这一划分在Linearizer.__init__中对应q_i/q_d独立/从属坐标与u_i/u_d独立/从属速度四个向量维度元组(l, m, n, o, s, k)会在构造时一次性推导出来见 linearize.py 的self._dims。一般地你可以任选哪些坐标/速度作为从属变量但实践中某些选择会在特定工作点引发数值奇异。如何系统地决定从属变量划分超出本文范围文档建议参考 Blajer (1994) 的坐标划分方法。线性化方程的两种输出形式sympy.physics.mechanics提供两种线性化 EOM 的输出形式。形式一M、A、B隐式形式默认在这种形式中力矩阵被线性化为两个独立矩阵 $A$ 和 $B$这是默认的线性化输出。方程为$$ M \begin{bmatrix} \delta \dot{q} \ \delta \dot{u} \ \delta \lambda \end{bmatrix} A \begin{bmatrix} \delta q_i \ \delta u_i \end{bmatrix} B \begin{bmatrix} \delta r \end{bmatrix} $$其中$$ \begin{aligned} M \in \mathbb{R}^{(nok) \times (nok)}\ A \in \mathbb{R}^{(nok) \times (n-lo-m)}\ B \in \mathbb{R}^{(nok) \times s} \end{aligned} $$注意 $q_i$ 和 $u_i$ 仅指独立坐标与速度而 $q$ 和 $u$ 同时包含独立与从属变量。这种形式保留了解线性方程组的自由度——你可以先符号化得到 M、A、B之后再代入数值工作点手工求解是大型符号系统的推荐做法。形式二A、B显式一阶状态空间形式这种形式把线性化 EOM 化为仅含独立坐标与速度的显式一阶形式常用于稳定性分析或控制理论$$ \begin{bmatrix} \delta \dot{q_i} \ \delta \delta \dot{u_i} \end{bmatrix} A \begin{bmatrix} \delta q_i \ \delta u_i \end{bmatrix} B \begin{bmatrix} \delta r \end{bmatrix} $$其中$$ \begin{aligned} A \in \mathbb{R}^{(n-lo-m) \times (n-lo-m)}\ B \in \mathbb{R}^{(n-lo-m) \times s} \end{aligned} $$使用该形式时在linearize方法中设置A_and_BTrue即可。对 Kane 方程线性化单摆示例下面用一个简单单摆系统演示完整流程。首先用KanesMethod建立系统并生成F_r与F_r^* from sympy import symbols, Matrix from sympy.physics.mechanics import * q1 dynamicsymbols(q1) # 摆的角度 u1 dynamicsymbols(u1) # 角速度 q1d dynamicsymbols(q1, 1) L, m, t, g symbols(L, m, t, g) # 构造世界坐标系 N ReferenceFrame(N) pN Point(N*) pN.set_vel(N, 0) # A.x 沿摆杆方向 A N.orientnew(A, axis, [q1, N.z]) A.set_ang_vel(N, u1*N.z) # 相对原点 N* 定位点 P P pN.locatenew(P, L*A.x) vel_P P.v2pt_theory(pN, N, A) pP Particle(pP, P, m) # 运动学微分方程 kde Matrix([q1d - u1]) # 作用在 P 点的力此处文档示例为水平力 R m*g*N.x # 用 Kane 法求解运动方程 KM KanesMethod(N, q_ind[q1], u_ind[u1], kd_eqskde) fr, frstar KM.kanes_equations([pP], [(P, R)])方式一直接使用 Linearizer 类通过to_linearizer类方法创建 linearizer 对象。该方法会把KanesMethod对象中的表示转换为上文描述的广义形式。由于独立/从属坐标与速度在创建KanesMethod时已经指定这里无需再次指定 linearizer KM.to_linearizer()然后调用Linearizer对象的linearize方法生成线性化 EOM M, A, B linearizer.linearize() M Matrix([ [1, 0], [0, -L**2*m]]) A Matrix([ [ 0, 1], [L*g*m*cos(q1(t)), 0]]) B Matrix(0, 0, [])也可以通过指定A_and_BTrue生成 A、B 形式 A, B linearizer.linearize(A_and_BTrue) A Matrix([ [ 0, 1], [-g*cos(q1(t))/L, 0]]) B Matrix(0, 0, [])工作点可以指定为字典或字典的可迭代对象。系统会在返回矩阵之前先在该指定点上求值 op_point {q1: 0, u1: 0} A_op, B_op linearizer.linearize(A_and_BTrue, op_pointop_point) A_op Matrix([ [ 0, 1], [-g/L, 0]])同样的效果也可以通过把msubs应用到未指定op_point时生成的矩阵上来实现 assert msubs(A, op_point) A_op有时返回的矩阵未必处于最简形式。你可以事后再做化简或者把Linearizer对象的simplify关键字参数设为True让它在内部执行化简。方式二使用 linearize 类方法KanesMethod类还提供了linearize方法作为便捷封装它内部调用to_linearizer执行线性化并返回结果。上文linearize方法中可用的所有关键字参数在这里同样可用 A, B, inp_vec KM.linearize(A_and_BTrue, op_pointop_point) A Matrix([ [ 0, 1], [-g/L, 0]])从源码看kane.py 中的KanesMethod.linearize实现为三行核心逻辑先self.to_linearizer(linear_solverlinear_solver)再linearizer.linearize(**kwargs)最后把输入向量linearizer.r追加到结果元组中见 kane.py。额外的输出inp_vec是一个包含所有未被广义坐标或速度向量收录的dynamicsymbols的向量。这些符号被视为系统的输入构成背景部分描述的 $r$ 向量。本例中没有输入所以向量为空 inp_vec Matrix(0, 0, [])从源码结构看r向量的构造发生在to_linearizer内部KanesMethod用find_dynamicsymbols在动力学驱动项f_3中搜索不属于q、qdot、u、udot、辅助速度及其导数的所有动态符号并用default_sort_key按正则序排序从而得到确定性的输入向量见 kane.py。注意早期 SymPy 版本中的KanesMethod线性化方法已被弃用将在未来版本移除。当前仓库中KanesMethod.linearize保留了new_method关键字参数用于向后兼容从源码看它已标注为 Deprecated, does nothing and will be removed见 kane.py新代码无需再显式传new_methodTrue。对 Lagrange 方程线性化Lagrange 方程的线性化流程与 Kane 方程基本相同。同样以单摆为例先把 A 坐标系与 P 点用q1d而非u1重新定义 # 用 q1d 而非 u1 重新定义 A 和 P A N.orientnew(A, axis, [q1, N.z]) A.set_ang_vel(N, q1d*N.z) P pN.locatenew(P, L*A.x) vel_P P.v2pt_theory(pN, N, A) pP Particle(pP, P, m) # 用 Lagrange 法求解运动方程 Lag Lagrangian(N, pP) LM LagrangesMethod(Lag, [q1], forcelist[(P, R)], frameN) lag_eqs LM.form_lagranges_equations()方式一直接使用 Linearizer 类LagrangesMethod对象同样可以用to_linearizer类方法生成Linearizer对象。与KanesMethod的唯一区别在于LagrangesMethod对象内部并没有预先指定独立/从属坐标与速度这些必须在调用to_linearizer时显式给出。本例没有从属坐标和速度如果有它们会包含在q_dep和qd_dep关键字参数中 linearizer LM.to_linearizer(q_ind[q1], qd_ind[q1d])此后流程与KanesMethod示例完全一致 A, B linearizer.linearize(A_and_BTrue, op_pointop_point) A Matrix([ [ 0, 1], [-g/L, 0]])方式二使用 linearize 类方法与KanesMethod类似LagrangesMethod类也提供linearize方法作为封装内部调用to_linearizer、执行线性化并返回结果。同样地独立/从属坐标与速度必须在调用时给出 A, B, inp_vec LM.linearize(q_ind[q1], qd_ind[q1d], A_and_BTrue, op_pointop_point) A Matrix([ [ 0, 1], [-g/L, 0]])从源码结构看lagrange.py 中的LagrangesMethod.to_linearizer会把 Lagrange 方程的各部分映射到广义形式f_0 u、f_1 -u运动学方程、f_2 _term1质量矩阵项、f_3 -(_term2 _term4)驱动力项、f_4 -_term3Lagrange 乘子项对应 $f_4$并把hol_coneqs与coneqs分别映射为 $f_c$ 与 $f_v$、对其求时间导数得到 $f_a$见 lagrange.py。该函数还会严格校验从属变量划分q_dep数量必须等于构型约束数、qd_dep数量必须等于速度约束数且q_ind q_dep必须恰好是q的一个划分否则抛出ValueError见 lagrange.py。另外带 Lagrange 乘子的 Lagrange 系统还可用solve_multipliers在指定工作点求解乘子值这在 test_linearize.py 的 nonminimal 单摆测试中被用于补全工作点字典属于约束系统线性化工作流的一部分。Linearizer 的内部机制置换矩阵、系数矩阵与线性求解器上面的 API 背后Linearizer做了几件重要的事理解它们有助于在大型系统上做出正确决策。延迟初始化_setup构造Linearizer时只保存广义方程向量与各变量向量置换矩阵、块矩阵、系数矩阵的计算被推迟到首次调用linearize时才执行_setup_done标志位以提高构造速度见 linearize.py。置换矩阵_form_permutation_matrices利用permutation_matrix构造Pq、Pu以及组合置换矩阵perm_mat满足 $[q_{ind}, u_{ind}]^T P\cdot[q, u]^T$。从属坐标/速度的排序约定是独立变量在前、从属变量在末尾P_qi/P_qd等子块按l、m切分见 linearize.py。系数矩阵 C_0、C_1、C_2_form_coefficient_matrices处理约束的影响。有构型约束$l 0$时 $C_0 (I - P_{qd}\cdot\mathrm{solver}(f_{c,jac},P_{qd}, f_{c,jac}))P_{qi}$无约束时 $C_0 I$。有运动约束$m 0$时类似地构造 $C_1$、$C_2$否则 $C_1 0$、$C_2 I$见 linearize.py。这一步正是正确处理约束的落点——约束导致从属变量与独立变量耦合一阶 Taylor 展开不能简单对右端求 Jacobian。块矩阵拼装_form_block_matrices对 $f_0, f_1, f_2, f_3, f_4$ 及各约束方程求 Jacobian按存在性条件拼出 $M$、$A$、$B$ 的各个子块如M_qq f_0.jacobian(qd)、A_qu -f_1.jacobian(u)等见 linearize.py。linearize方法再按注释中的块结构把这些子块row_join成完整矩阵。可插拔线性求解器linear_solver参数既可以是MatrixBase.solve接受的方法名字符串默认LU对应A.LUsolve(b)也可以是形如x f(A, b)的可调用对象。_parse_linear_solver负责解析见 functions.py。文档特别提示LUsolve计算快但常常遇到零除而产生nan结果——这正是排障章节第二个问题的底层原因之一。linearize方法的关键字参数汇总参数默认值说明op_pointNone工作点字典或字典的可迭代对象值可为数值或符号符号被替换得越多运行越快A_and_BFalseFalse返回(M, A, B)True返回(A, B)simplifyFalse是否在返回前对矩阵做符号化简对大表达式可能非常耗时潜在问题与排障虽然Linearizer类理论上能线性化所有系统实践中仍会遇到一些问题。问题一符号化线性化且 A_and_BTrue 时很慢最可能的原因是符号化求解大型线性方程组本身是昂贵操作。指定工作点可以缩小表达式规模、加速求解。如果确实需要纯符号解例如之后要批量代入许多工作点可以绕过A_and_BTrue的内部符号求解先用A_and_BFalse取回 M、A、B代入工作点后再手工求解 M, A, B linearizer.linearize() M_op msubs(M, op_point) A_op msubs(A, op_point) perm_mat linearizer.perm_mat A_lin perm_mat.T * M_op.LUsolve(A_op) A_lin Matrix([ [ 0, 1], [-g/L, 0]])这与源码注释给出的公式一致$A P^T M^{-1} A$、$B P^T M^{-1} B$其中 $P \texttt{Linearizer.perm_mat}$见 linearize.py。求解前A、M中符号越少求解越快对大型表达式把转换到 A、B 形式推迟到大部分符号已代入数值之后收益明显。问题二线性化结果出现 nan、zoo 或 oo有两种可能的成因。成因一优先检查从属坐标选择在特定工作点导致奇异。某些坐标划分在特定工作点会产生奇异。系统性地避开它需要坐标划分理论超出本文范围参见 Blajer 1994。成因二代入工作点前矩阵未充分化简。一个典型的演示 from sympy import sin, tan expr sin(q1)/tan(q1) op_point {q1: 0} expr.subs(op_point) nan如果先化简再代入就能得到正确值 expr.simplify().subs(op_point) 1目前尚未找到完全自动避免该问题的通用办法。对规模合理的表达式用msubs并设smartTrue会启用一种规避此类条件的算法但对大型表达式它极其耗时 msubs(expr, op_point, smartTrue) 1从源码结构看msubs是physics.mechanics专用的自研替换函数见 functions.py它只遍历表达式树一次跳过Derivative内部项并且支持传入多个替换字典smartTrue时会额外检查直接代入会产生nan、但化简后有效的条件。验证与进阶测试套件中的非最小坐标示例上述单摆示例虽简单但没有包含从属坐标或从属速度。仓库测试 test_linearize.py 提供了可直接运行的更复杂示例nonminimal 单摆Kane 法用直角坐标q1, q2摆锤的 N.x、N.y 坐标与速度u1, u2建模显式给出构型约束f_c |P.pos_from(pN)| - L与速度约束f_v P.vel(N).express(A).dot(A.x)通过q_dependent[q1], u_dependent[u1]指定从属变量在直立工作点线性化后得到A [[0, 1], [-9.8/L, 0]]见 test_linearize.py。该测试同时验证了linear_solverGJ与自定义 callable 求解器三种方式的等价性。nonminimal 单摆Lagrange 法用hol_coneqsf_c建立约束先用solve_multipliers求出 Lagrange 乘子的工作点值并入op_point再线性化验证见 test_linearize.py。滚动圆盘Kane 法含加速度约束6 个广义坐标、6 个广义速度的非完整约束系统在直立滚动名义点线性化后验证 A 矩阵并检查临界速度处特征值全为零见 test_linearize.py。这个用例展示了op_point接收多字典列表[q_op, u_op, qd_op, ud_op]的用法。此外文档还指引读者参考 Tutorials 页中 Nonminimal Coordinates Pendulum 一节其中用 Kane 与 Lagrange 两种方法对同一单摆分别带从属坐标做了线性化可作为本指南示例之外的延伸阅读。小结SymPy 的线性化路径可以概括为一条清晰的调用链KanesMethod/LagrangesMethod生成 EOM →to_linearizer强制转换为 $f_c/f_v/f_a/f_0\sim f_4$ 广义形式自动完成约束感知→Linearizer.linearize做一阶 Taylor 展开并按块矩阵拼装。无约束时它退化为右端对 $q$、$u$ 的 Jacobian有约束时通过置换矩阵与 $C_0, C_1, C_2$ 系数矩阵正确处理从属变量。实践要点大型符号系统优先用M/A/B形式加op_point数值代入再手工求解出现nan时先怀疑坐标划分其次用msubs(..., smartTrue)或化简后再代入。相关实现与测试分别位于 linearize.py、kane.py、lagrange.py 与 test_linearize.py可按需深入。【免费下载链接】sympyA computer algebra system written in pure Python项目地址: https://gitcode.com/GitHub_Trending/sy/sympy创作声明：本文部分内容由AI辅助生成（AIGC），仅供参考