SymPy 多体系统符号线性化指南：Linearizer 类的原理、用法与约束处理

SymPy 多体系统符号线性化指南Linearizer 类的原理、用法与约束处理【免费下载链接】sympyA computer algebra system written in pure Python项目地址: https://gitcode.com/GitHub_Trending/sy/sympy本文围绕 SymPyphysics.mechanics模块中的Linearizer类docstring 页面见 linearize.rst完整实现位于 linearize.py展开系统讲解如何在多体动力学系统建模后求解其线性化状态空间方程重点覆盖一般系统形式general form的组织方式、配置约束与速度约束下独立/相关坐标的处理、操作点operating point的设置技巧以及线性求解器的选择。读完本文你将能够独立完成从KanesMethod/LagrangesMethod建模到A, B状态矩阵提取的完整符号线性化流程。一、Linearizer 是什么约束多体系统的线性化求解器在控制设计与稳定性分析中我们常常需要将非线性的运动学方程与动力学方程在某平衡点附近线性化。Linearizer正是 SymPy 中负责这一任务的类它持有动态系统的一般模型形式用于计算系统的线性化形式并且能够正确处理由约束引入的相关坐标dependent coordinates与相关速度dependent speeds。其核心算法与符号约定完整描述于 Peterson、Gede 与 Hubbard 的论文Symbolic linearization of equations of motion of constrained multibody systemsMultibody Syst Dyn, 2015中docstring 中给出了该文献引用。Linearizer的构造函数位于 linearize.py签名如下Linearizer(f_0, f_1, f_2, f_3, f_4, f_c, f_v, f_a, q, u, q_iNone, q_dNone, u_iNone, u_dNone, rNone, lamsNone, linear_solverLU)在实际使用中几乎不需要手工构造这些f矩阵——推荐通过KanesMethod.to_linearizer()见 kane.py或LagrangesMethod.to_linearizer()见 lagrange.py自动从运动方程中组装出Linearizer实例。直接构造该类的场景仅限于用户希望完全自定义系统形式、或需要复用同一个Linearizer在多个操作点反复线性化以提升效率时。二、一般系统形式General Form与各矩阵的含义Linearizer将系统运动方程组织为如下一般形式这也是f_0~f_4八个矩阵的语义来源f_c(q, t) 0 配置完整约束 f_v(q, u, t) 0 速度非完整约束 f_a(q, u, u, t) 0 加速度约束f_v 的时间导数 f_0(q, q, u, t) 0 运动学微分方程中不含 u 的部分 f_1(q, q, u, t) 0 运动学微分方程中与 u 相关的部分 f_2(q, u, u, r, t) 0 动力学方程中与 u 相关的部分 f_3(q, u, u, r, t) 0 动力学方程中其余部分含主动力 f_4(q, u, u, lams, t) 0 与拉格朗日乘子 lams 相关的部分这些矩阵在 linearize.py 中被统一转换为Matrix存储并公开为同名属性f_0~f_4、f_c、f_v、f_a。对于不存在的项传入空数组或空 Matrix即可。从KanesMethod.to_linearizer()的实现kane.py可以看出这些矩阵的组装逻辑f_c取自身份self._f_h完整约束f_v由非完整约束与速度约束系数矩阵组装f_a为f_v对时间求导f_0与f_1由运动学微分方程self._f_k分别对u0与q0做msub拆分得到f_2与f_3由广义惯性力frstar与广义主动力fr拆分得到f_4对 Kane 法而言为零矩阵而对拉格朗日法lagrange.py而言f_4 -self._term3承载与拉格朗日乘子相关的项。变量向量与维度参数构造时传入的变量向量也被保存为属性属性含义q广义坐标向量u广义速度向量r输入变量forcing inputs向量lams拉格朗日乘子向量q_i/u_i独立广义坐标 / 独立广义速度q_d/u_d相关dependent广义坐标 / 相关广义速度perm_mat置换矩阵满足[q_ind, u_ind]^T perm_mat * [q, u]^T其中r的自动识别规则是运动方程中所有不属于q、u、q、u的动态符号dynamicsymbols按规范序排序后构成输入向量见 kane.py。同时源码会对r中同时出现变量及其导数的情形抛出ValueError避免对指定量的导数做线性化。在 linearize.py 中Linearizer还会根据这些向量的长度推导出六个维度参数并以命名元组_dims (l, m, n, o, s, k)保存l配置约束的个数m速度约束的个数n广义坐标个数o广义速度个数s输入变量个数k拉格朗日乘子个数。这些维度参数贯穿整个线性化过程决定了后面各分块矩阵是否存在。三、约束处理置换矩阵与系数矩阵 C_0、C_1、C_2约束系统线性化的核心难点在于区分独立/相关坐标与速度。Linearizer通过**延迟初始化lazy setup**机制处理这一问题昂贵的矩阵构造被放在_setup()linearize.py中仅在第一次调用linearize()时执行一次从而加快Linearizer对象的创建速度。置换矩阵 Pq、Pu_form_permutation_matrices()linearize.py调用模块级辅助函数permutation_matrix(orig_vec, per_vec)linearize.py构造置换矩阵。该辅助函数要求两个向量长度相同且包含完全相同的符号否则抛出ValueError返回的置换矩阵P满足orig_vec P * per_vec。由此得到_Pq将q重排为[q_i, q_d]的置换矩阵进一步切分为_Pqi对应独立坐标列与_Pqd对应相关坐标列仅当l 0时存在_Pu同理处理u与[u_i, u_d]切分为_Pui与_Pudperm_mat由_Pqi、_Pui组合而成的整体置换矩阵用于从全状态[q, u]中提取独立状态[q_i, u_i]。系数矩阵 C_0、C_1、C_2_form_coefficient_matrices()linearize.py通过求解约束雅可比矩阵的线性系统构造三个关键系数矩阵C_0当存在配置约束l 0时C_0 (I_n - P_qd * (f_c_q * P_qd)^{-1} * f_c_q) * P_qi其中f_c_q f_c.jacobian(q)若无配置约束则退化为I_nn 阶单位阵C_1当存在速度约束m 0时C_1 -P_ud * (f_v_u * P_ud)^{-1} * f_v_q其中f_v_u f_v.jacobian(u)、f_v_q f_v.jacobian(q)否则为零矩阵zeros(o, n)C_2有速度约束时为(I_o - P_ud * (f_v_u * P_ud)^{-1} * f_v_u) * P_ui无速度约束时为I_o。这三个矩阵的物理含义是把相关坐标/速度对独立坐标/速度的依赖关系由约束隐式定义显式地吸收进线性化系数从而将约束系统的状态约化到独立子空间。四、linearize() 方法从隐式到显式的两种输出形式linearize()是核心入口方法linearize.py签名如下linearize(op_pointNone, A_and_BFalse, simplifyFalse)参数说明参数类型/默认值说明op_pointdict或Iterable[dict]默认None操作点条件覆盖广义坐标、广义速度及其对时间导数的全部或子集。None表示操作点保持为任意符号。提前代入的符号越多运行越快A_and_Bbool默认FalseFalse时返回(M, A, B)隐式形式True时返回(A, B)显式状态空间形式simplifybool默认False返回前是否对结果调用simplify()。表达式很大时可能非常耗时返回矩阵与系统方程当A_and_BFalse默认时返回矩阵M, A, B对应隐式形式[M] * [q, u]^T [A] * [q_ind, u_ind]^T [B] * r当A_and_BTrue时返回矩阵A, B对应显式状态空间形式[q_ind, u_ind]^T [A] * [q_ind, u_ind]^T [B] * r从 linearize.py 的实现可以看到M、A、B均按分块矩阵组装质量矩阵 M三行分块[M_qq; M_uqc; M_uqd]、[0; M_uuc; M_uud]、[0; 0; M_uld]即M |M_qq 0_nxo 0_nxk | |M_uqc M_uuc 0_mxk | |M_uqd M_uud M_uld |状态系数矩阵 A由A_qq, A_qu, A_uqc, A_uuc, A_uqd, A_uud与C_0, C_1, C_2组合即A |(A_qq A_qu*C_1)*C_0 A_qu*C_2 | |(A_uqc A_uuc*C_1)*C_0 A_uuc*C_2 | |(A_uqd A_uud*C_1)*C_0 A_uud*C_2 |输入矩阵 BB |0_(nm)xs; B_u|仅当s ! 0时存在否则返回空矩阵。所有矩阵在组装后都会通过msub将op_point字典代入源码使用msub见 linearize.py这也是提前代入符号可加速的原因。分块矩阵_M_qq、_A_qq、_M_uqc等的构造位于_form_block_matrices()linearize.py它们分别由f_0~f_4对q、u、q、u、lams、r求雅可比得到。关于 A_and_BTrue 的计算成本与替代方案docstring 特别提醒若系统含有大量符号参数A_and_BTrue需要求解M*x A、M*x B两组符号线性系统linearize.py计算量很大。此时更推荐使用默认的A_and_BFalse先拿到M, A, B事后将数值代入再按下式恢复显式状态空间形式A P.T * M.LUsolve(A) B P.T * M.LUsolve(B)其中P Linearizer.perm_mat。这正是 docstring Notes 中给出的标准做法。五、linear_solver 参数字符串与可调用对象Linearizer与KanesMethod/LagrangesMethod的线性化接口都接受linear_solver参数用于求解线性化过程中反复出现的符号线性系统A*x b如f_c_q * P_qd的求逆。该参数由辅助函数_parse_linear_solver()functions.py解析def _parse_linear_solver(linear_solver): if callable(linear_solver): return linear_solver return lambda A, b: Matrix.solve(A, b, methodlinear_solver)字符串必须是MatrixBase.solve()支持的合法方法名底层转为A.solve(b, method...)可调用对象需符合x f(A, b)接口直接原样返回使用默认值LU对应 SymPy 的A.LUsolve(b)。docstring 明确指出LUsolve()计算快但经常因除零导致nan结果——这是符号线性化实践中常见的坑。因此当默认LU出现除零/nan时可以尝试切换为GJGauss-Jordan 消元法或传入自定义可调用对象。测试 test_linearize.py 验证了linear_solverGJ与默认LU得到完全一致的线性化结果test_functions.py 则直接断言了_parse_linear_solver对可调用对象与字符串两种形态的解析行为。测试中还注明 symengine 后端不支持methodGJ这是选用求解器时需要注意的兼容性限制。六、从 KanesMethod 与 LagrangesMethod 出发的两种工作流Linearizer通常不直接实例化而是作为KanesMethod与LagrangesMethod高层接口的底层引擎。两条路径的推荐流程如下。6.1 高层便捷接口直接调用 linearize()KanesMethod.linearize()kane.py在(M, A, B)或(A, B)之外额外返回第三个元素r输入向量即A, B, r或M, A, B, rr为方程中不属于q, u, q, u的动态符号按规范序排序LagrangesMethod.linearize(q_ind, qd_ind, q_dep, qd_dep, ...)lagrange.py则需要显式指定独立/相关的坐标与速度划分。两者的linear_solver参数语义与Linearizer完全一致其余关键字参数op_point、A_and_B、simplify会原样透传给Linearizer.linearize()。6.2 可复用对象to_linearizer()如果需要在多个操作点反复线性化两个类都提供了to_linearizer()它一次性组装Linearizer对象后续每次linearizer.linearize(op_point...)只做代入与求解比重建整个方程再线性化高效得多。docstring 明确建议这种方式见 kane.py、lagrange.py。七、完整实战示例非最小实现摆的线性化下面以测试 test_linearize.py 中的非最小实现non-minimal摆为例演示包含配置约束与速度约束的完整线性化流程。该模型用平面坐标(q1, q2)描述摆锤位置摆长固定构成约束摆沿杆方向速度恒为零构成速度约束。from sympy import symbols, Matrix, atan, solve from sympy.physics.mechanics import (dynamicsymbols, ReferenceFrame, Point, Particle, KanesMethod) q1, q2 dynamicsymbols(q1:3) # 平面坐标 q1d, q2d dynamicsymbols(q1:3, level1) u1, u2 dynamicsymbols(u1:3) u1d, u2d dynamicsymbols(u1:3, level1) L, m, t symbols(L, m, t) g 9.8 N ReferenceFrame(N) pN Point(N*) pN.set_vel(N, 0) # A.x 沿摆方向 theta1 atan(q2/q1) A N.orientnew(A, axis, [theta1, N.z]) P pN.locatenew(P1, q1*N.x q2*N.y) pP Particle(pP, P, m) # 运动学微分方程 kde Matrix([q1d - u1, q2d - u2]) dq_dict solve(kde, [q1d, q2d]) P.set_vel(N, P.pos_from(pN).dt(N).subs(dq_dict)) # 配置约束摆长恒为 L f_c Matrix([P.pos_from(pN).magnitude() - L]) # 速度约束沿杆方向速度恒为零 f_v Matrix([P.vel(N).express(A).dot(A.x)]) # 加速度约束速度约束的时间导数 f_a f_v.diff(t) # 重力载荷 R m*g*N.x # 建立 Kane 方程q1 相关、u1 相关q2/u2 独立 KM KanesMethod(N, q_ind[q2], u_ind[u2], q_dependent[q1], u_dependent[u1], configuration_constraintsf_c, velocity_constraintsf_v, acceleration_constraintsf_a, kd_eqskde) (fr, frstar) KM.kanes_equations([pP], [(P, R)]) # 定义操作点摆竖直向下、静止 q_op {q1: L, q2: 0} u_op {u1: 0, u2: 0} ud_op {u1d: 0, u2d: 0} # 线性化得到显式状态空间矩阵 A, B, inp_vec KM.linearize(op_point[q_op, u_op, ud_op], A_and_BTrue, simplifyTrue) # 期望结果对应线性化摆方程 [0 1; -g/L 0] assert A.expand() Matrix([[0, 1], [-9.8/L, 0]]) assert B Matrix([])要点解读操作点可以是字典的列表op_point[q_op, u_op, ud_op]被合并为一个字典后整体代入这正是 linearize.py 中对Iterable类型的处理逻辑约束处理是自动的q1/u1被声明为相关坐标/速度后C_0、C_1、C_2会自动吸收约束最终得到 2 阶状态矩阵A状态为独立量q2, u2且稳态摆的线性化矩阵与经典结果[[0, 1], [-g/L, 0]]完全一致inp_vecKanesMethod.linearize额外返回的输入向量r本例无外部时变输入故B为空矩阵。八、复杂案例滚动圆盘含六维约束与拉格朗日路线测试文件中还提供了更复杂的验证案例——受约束滚动圆盘rolling disc。在test_linearize_rolling_disc_kanetest_linearize.py中圆盘有 6 个广义坐标q1:q6其中q6为相关坐标满足配置约束f_c [q6 - dot(CO.pos_from(P), N.z)]圆盘圆心高度与接触几何一致6 个广义速度中u4, u5, u6为相关速度满足速度约束f_v [dot(P.vel(N), uv) for uv in C]接触点速度为零即纯滚动无滑动条件操作点由q_op, u_op, qd_op, ud_op四个字典组成代入后线性化最终验证了 8 阶矩阵A在竖直稳态标称点处的取值并验证了在临界速度q3d 1/sqrt(3)下所有特征值为 0 的稳定性结论。同一物理系统也可走拉格朗日路线test_linearize_rolling_disc_lagrangetest_linearize.py先由Lagrangian构造LagrangesMethod再以l.linearize(q_indq, qd_indqd, op_pointop_point, A_and_BTrue)得到 6 阶状态矩阵与 Kane 路线结果相互印证。这体现了Linearizer作为统一底层引擎对两种建模方法Kane 法与 Lagrange 法的通用支撑。九、工程实践建议与常见陷阱优先用隐式形式A_and_BFalse符号参数较多时A_and_BTrue需要求解大规模符号线性系统开销巨大先取M, A, B、代入数值后再用A perm_mat.T * M.LUsolve(A)恢复显式形式是更稳健的路径。警惕LU的除零问题docstring 明确警告LUsolve()可能因除零产生nan。遇到时切换到GJ或自定义求解器如lambda A, b: A.LUsolve(b)但需注意 symengine 后端不支持methodGJ。操作点代入越充分越快op_point中给出的符号替换越多后续符号表达式越小、求解越快操作点必须满足运动方程即q_op, u_op, qd_op, ud_op须为实际平衡/稳态点。拉格朗日法需显式划分独立/相关量LagrangesMethod.linearize(q_ind, qd_ind, q_dep, qd_dep, ...)要求独立与相关划分与约束个数严格匹配len(q_dep) len(f_c)、len(u_dep) len(f_v)且q_ind q_dep必须恰好等于全部坐标否则抛出ValueError见 lagrange.py。输入向量 r 中禁止出现其导数to_linearizer()会检查r中是否存在变量的时间导数若同时出现则抛出ValueError因为线性化强制项无法处理这类耦合。通过以上内容你可以基于 Linearizer 源码、KanesMethod 与 LagrangesMethod 的接口以及 test_linearize.py 中的完整验证案例将任意带约束的多体系统高效地符号线性化为后续的控制设计、特征值分析与稳定性判断奠定基础。【免费下载链接】sympyA computer algebra system written in pure Python项目地址: https://gitcode.com/GitHub_Trending/sy/sympy创作声明：本文部分内容由AI辅助生成（AIGC），仅供参考