尧图网站设计 尧图网站设计YAOTU DESIGN
ARTICLE DETAIL

资讯详情

深耕网站设计与一线实操的经验洞察。

Sympy SymbolicSystem:统一表达多体系统运动方程的三种标准形式

Sympy SymbolicSystem:统一表达多体系统运动方程的三种标准形式 Sympy SymbolicSystem统一表达多体系统运动方程的三种标准形式【免费下载链接】sympyA computer algebra system written in pure Python项目地址: https://gitcode.com/GitHub_Trending/sy/sympy本文基于 Sympy 官方解释文档《Symbolic Systems in Physics/Mechanics》完整讲解sympy.physics.mechanics中的SymbolicSystem类它以统一的数据格式承载多体动力系统的运动方程支持显式与隐式三种方程形态并可选携带系统刚体、载荷、代数约束行号等元信息。读完本文你将能够把任意多体系统含约束的 DAE 系统的运动方程手工输入到SymbolicSystem中正确区分三种方程形式与参数约定并通过属性与compute_explicit_form()在形式之间安全转换为后续对接数值求解器ODE/DAE 求解代码做好准备。一、SymbolicSystem 的定位与设计目标SymbolicSystem是physics/mechanics模块中用于存放多体动力系统全部相关信息的容器类。按其文档说明源文件它的核心定位是最基本形态包含系统的运动方程equations of motion, EOM可选扩展信息系统所受载荷loads、组成系统的刚体/质点bodies以及用户认为重要的任何附加方程设计目标为运动方程提供统一的输出格式使数值分析代码可以围绕这一格式进行设计。也就是说SymbolicSystem并不负责“从模型推导方程”而是负责把已经得到的方程无论是解析推导还是其他工具生成规范地组织起来供数值积分、线性化等下游环节消费。需要与同模块中另一个类System定义区分开System是面向建模的高层类通过add_bodies、add_joints、add_loads、add_kdes等接口搭建模型再由form_eoms(eom_methodKanesMethod, **kwargs)实现调用KanesMethod或LagrangesMethod后端生成运动方程而SymbolicSystem是面向方程本身的数据容器两者都通过 模块导出 对外暴露可以from sympy.physics.mechanics import SymbolicSystem或import sympy.physics.mechanics.system as system两种方式引入。二、三种运动方程标准形式SymbolicSystem支持三种等价的输入形式这是理解整个类的关键。设x为状态量如[q, u]t为时间r为指定的外部输入exogenous inputsp为常数q为广义坐标u为广义速度形式结构方程说明[1] 显式合并形式运动学与动力学合并、显式x F_1(x, t, r, p)直接给出全部状态的导数[2] 隐式合并形式运动学与动力学合并、隐式M_2(x, p) x F_2(x, t, r, p)左侧为质量矩阵[3] 隐式分离形式运动学与动力学分开、隐式M_3(q, p) u F_3(q, u, t, r, p)q G(q, u, t, r, p)动力学隐式、运动学显式其中各符号的含义为F_1合并方程显式形式的右端F_2合并方程隐式形式的右端F_3动力学方程隐式形式的右端M_2合并方程隐式形式的质量矩阵M_3动力学方程隐式形式的质量矩阵即广义惯量阵G运动学微分方程的右端。从源码看__init__依据传入参数自动判定形式形式判定逻辑若提供了coordinate_derivatives即G判为形式 [3]此时right_hand_side被解释为F_3mass_matrix被解释为M_3否则若提供了mass_matrix判为形式 [2]right_hand_side被解释为F_2mass_matrix为M_2否则判为形式 [1]right_hand_side即F_1。这一“由参数组合决定解释方式”的机制是后续所有属性可用性与报错行为的根源。三、完整实例单摆的笛卡尔坐标模型手工输入下面的示例与官方文档 symsystem.rst 保持一致以笛卡尔坐标描述单摆质量点位置作为广义坐标而非通常的最小坐标形式把运动方程手工输入SymbolicSystem。该模型与 lin_pend_nonmin_example 教程中的非最小坐标单摆等价——教程使用q1、q2作为质量点的水平/竖直坐标而本文示例将其替换为x、y且参考系相对教程旋转了 90 度因此重力载荷沿N.x方向。3.1 导入与符号初始化from sympy import atan, symbols, Matrix from sympy.physics.mechanics import (dynamicsymbols, ReferenceFrame, Particle, Point) import sympy.physics.mechanics.system as system from sympy.physics.vector import init_vprinting init_vprinting(pretty_printFalse) # 动态符号时间的函数位置 x, y广义速度 u, v约束乘子 lam x, y, u, v, lam dynamicsymbols(x y u v lambda) # 常数符号质量、摆长、重力加速度 m, l, g symbols(m l g)状态向量为(x, y, u, v, lam)共 5 个分量其中lam是用于约束x**2 y**2 l**2的乘子——这正是方程呈微分代数方程DAE形式的原因。3.2 用三种形式写出同一套运动方程先定义形式 [3] 的动力学方程M_3 u F_3与形式 [2] 的合并隐式方程M_2 x F_2# 形式 [3]3x3 质量矩阵 右端覆盖 [u, v, lam] dyn_implicit_mat Matrix([[1, 0, -x/m], [0, 1, -y/m], [0, 0, l**2/m]]) dyn_implicit_rhs Matrix([0, 0, u**2 v**2 - g*y]) # 形式 [2]5x5 块矩阵 右端前两行即运动学方程 x u, y v comb_implicit_mat Matrix([[1, 0, 0, 0, 0], [0, 1, 0, 0, 0], [0, 0, 1, 0, -x/m], [0, 0, 0, 1, -y/m], [0, 0, 0, 0, l**2/m]]) comb_implicit_rhs Matrix([u, v, 0, 0, u**2 v**2 - g*y]) # 形式 [3] 的运动学右端 G kin_explicit_rhs Matrix([u, v]) # 形式 [1]对隐式系统做 LU 分解求解得到显式右端 comb_explicit_rhs comb_implicit_mat.LUsolve(comb_implicit_rhs)注意形式 [2] 的comb_implicit_mat具有块对角结构左上块前两行是运动学方程的隐式矩阵单位阵右下块后三行就是M_3。源码在从形式 [3] 自动构造comb_implicit_mat属性时正是按此结构做零块拼接见下文 5.2 节。3.3 建立参考系、质点与载荷虽然方程已手工给出仍可建立刚体/质点对象以便将bodies与loads一并传入容器theta atan(x/y) omega dynamicsymbols(omega) N ReferenceFrame(N) A N.orientnew(A, Axis, [theta, N.z]) A.set_ang_vel(N, omega * N.z) O Point(O) O.set_vel(N, 0) P O.locatenew(P, l * A.x) P.v2pt_theory(O, N, A) # 得到速度 l*omega*A.y Pa Particle(Pa, P, m) bodies [Pa] loads [(P, g * m * N.x)]载荷的书写约定为力用(作用点, 力矢量)的元组力矩用(受作用的参考系, 力矩矢量)的元组见 类文档字符串。3.4 标记代数约束行alg_con本示例的运动方程是DAE形式DAE 求解器需要知道哪些行是代数方程而非微分方程。这一信息通过alg_con以行索引列表传入SymbolicSystem有一个重要的索引约定传入时行索引对应你输入的矩阵——形式 [3] 下对应dyn_implicit_mat3 行本例为alg_con [2]访问时alg_con属性始终对应合并了运动学与动力学的完整方程5 行本例为alg_con_full [4]。alg_con [2] # 形式 [3] 下的输入索引 alg_con_full [4] # 合并形式下的索引源码在__init__中会自动完成这一偏移若同时提供了coordinate_derivatives且给出alg_con则对每个索引加上运动学方程行数偏移逻辑因此形式 [3] 传入[2]后symsystem.alg_con属性返回的就是[4]。3.5 状态向量与坐标/速度索引states (x, y, u, v, lam) coord_idxs (0, 1) # 前两个状态分量是广义坐标 speed_idxs (2, 3) # 第 3、4 个分量是广义速度当第一个参数coord_states传入的是完整状态而非单独的坐标时coord_idxs/speed_idxs告诉SymbolicSystem哪些分量是坐标、哪些是速度若不传入coordinates与speeds属性将无法访问并抛出AttributeError。3.6 创建三个等价的 SymbolicSystem 实例三种形式分别构造参数与形式的对应关系是本节重点# 形式 [1]只给显式右端 symsystem1 system.SymbolicSystem(states, comb_explicit_rhs, alg_conalg_con_full, bodiesbodies, loadsloads) # 形式 [2]右端 合并质量矩阵 symsystem2 system.SymbolicSystem(states, comb_implicit_rhs, mass_matrixcomb_implicit_mat, alg_conalg_con_full, coord_idxscoord_idxs) # 形式 [3]动力学右端 动力学质量矩阵 运动学右端 symsystem3 system.SymbolicSystem(states, dyn_implicit_rhs, mass_matrixdyn_implicit_mat, coordinate_derivativeskin_explicit_rhs, alg_conalg_con, coord_idxscoord_idxs, speed_idxsspeed_idxs)各实例属性的实际输出 symsystem1.states Matrix([ [x], [y], [u], [v], [lambda]]) symsystem2.coordinates Matrix([ [x], [y]]) symsystem3.speeds Matrix([ [u], [v]]) symsystem1.comb_explicit_rhs Matrix([ [u], [v], [(-g*y u**2 v**2)*x/l**2], [(-g*y u**2 v**2)*y/l**2], [m*(-g*y u**2 v**2)/l**2]]) symsystem2.comb_implicit_rhs Matrix([ [u], [v], [0], [0], [-g*y u**2 v**2]]) symsystem2.comb_implicit_mat Matrix([ [1, 0, 0, 0, 0], [0, 1, 0, 0, 0], [0, 0, 1, 0, -x/m], [0, 0, 0, 1, -y/m], [0, 0, 0, 0, l**2/m]]) symsystem3.dyn_implicit_rhs Matrix([ [0], [0], [-g*y u**2 v**2]]) symsystem3.dyn_implicit_mat Matrix([ [1, 0, -x/m], [0, 1, -y/m], [0, 0, l**2/m]]) symsystem3.kin_explicit_rhs Matrix([ [u], [v]]) symsystem1.alg_con [4] symsystem1.bodies (Pa,) symsystem1.loads ((P, g*m*N.x),)几个值得注意的行为symsystem3.alg_con同样返回[4]——尽管构造时传入的是[2]源码已自动完成偏移bodies与loads在构造时被统一转换为tuple存储转换代码所以symsystem1.bodies输出为(Pa,)与coord_idxs同理bodies/loads只有在初始化时指定过属性才可访问未指定则访问抛AttributeError。四、最小示例单自由度形式 [3] 输入类文档字符串示例部分给出了一个更简洁的用法单摆用角度theta作广义坐标、omega作广义速度直接按形式 [3] 输入坐标与速度通过位置参数分开传入此时第一个参数是坐标speeds作为第三个位置参数from sympy import Matrix, sin, symbols from sympy.physics.mechanics import dynamicsymbols, SymbolicSystem l, m, g symbols(l m g) theta, omega dynamicsymbols(theta omega) kin_explicit_rhs Matrix([omega]) # G dyn_implicit_mat Matrix([l**2 * m]) # M_3 dyn_implicit_rhs Matrix([-g * l * m * sin(theta)]) # F_3 symsystem SymbolicSystem([theta], dyn_implicit_rhs, [omega], dyn_implicit_mat)对比 3.6 节可知两种传参风格的差异传参风格第一个参数speeds适用场景坐标/速度分开广义坐标集合第三个位置参数speeds方程按最小坐标编写如本例完整状态全部状态[q, u, ...]不传状态中混有乘子等非速度量需coord_idxs/speed_idxs区分源码对两种风格的分支处理见 初始化逻辑传入speeds时states由coordinates.col_join(speeds)拼成不传时第一个参数整体作为states。五、源码级实现解析5.1 构造函数与形式判定SymbolicSystem的完整构造签名为源码def __init__(self, coord_states, right_hand_side, speedsNone, mass_matrixNone, coordinate_derivativesNone, alg_conNone, output_eqns{}, coord_idxsNone, speed_idxsNone, bodiesNone, loadsNone):参数要点coord_states有序可迭代时间的函数集合是否表示“坐标”取决于speeds是否提供right_hand_sideMatrix其语义F_1/F_2/F_3由mass_matrix、coordinate_derivatives是否传入决定speeds提供后第一参数被解释为广义坐标mass_matrix形式 [2]/[3] 的质量矩阵coordinate_derivatives提供即宣告形式 [3]内容是运动学方程q G的右端alg_con代数约束行索引形式 [3] 下按mass_matrix/right_hand_side的行编号解释内部自动偏移为合并形式的行号output_eqns字典键为输出方程名、值为符号表达式用于追踪额外要跟踪的输出量单摆示例中可为{PE: m*g*(ly)}见 测试coord_idxs/speed_idxs当第一参数是完整状态时指定坐标/速度对应的索引bodies/loads刚体对象与载荷力为(point, force)、力矩为(frame, torque)的可迭代对象内部转为 tuple。5.2 属性的自动推导与形式互转各矩阵属性并非全部需要用户提供源码按以下规则推导comb_implicit_mat/comb_implicit_rhs若以形式 [3] 输入访问这两个属性时源码会用eye(n_kin)、zeros块自动拼接出合并矩阵/右端实现与 3.2 节手工写的块矩阵结构完全一致comb_explicit_rhs不能直接访问必须先调用compute_explicit_form()计算实现。该方法内部即对隐式系统做LUsolve形式 [3] 下先解M_3 u F_3再与kin_explicit_rhs纵向拼接形式 [2] 下直接解M_2 x F_2。源码注释特别提示该计算“potentially take awhile”因此被设计为显式调用而非惰性自动计算states始终可用coordinates/speeds仅在初始化时能确定对应分量时可用只读约束所有属性均为只读 property尝试赋值如symsystem.bodies 42会抛AttributeError见 属性只读测试。5.3 未指定信息的错误行为测试test_not_specified_errors源码系统覆盖了越界访问的行为值得在集成时记住形式 [1] 的实例访问comb_implicit_mat、dyn_implicit_rhs、kin_explicit_rhs等均未定义抛AttributeError且因为comb_explicit_rhs已存在compute_explicit_form()也会报错守卫逻辑形式 [2] 的实例访问dyn_implicit_mat、kin_explicit_rhs报错——源码不会从合并矩阵中反拆出动力学块仅传状态而未传coord_idxs/speed_idxs时访问coordinates/speeds报错未传bodies/loads时访问对应属性报错。5.4 动态符号与常数符号的自动提取SymbolicSystem还提供两个辅助方法便于在对接数值求解器时确定变量清单dynamic_symbols()实现扫描运动方程表达式显式右端或隐式矩阵与右端的全部元素中的Dynamicsymbol再并入状态量返回时间依赖符号的 tupleconstant_symbols()实现取同样的表达式集合的free_symbols剔除时间符号t返回常数符号的 tuple。对本文示例两者分别得到{x, y, u, v, lam}与{m, l, g}与 测试断言 一致。六、API 速查表构造函数参数参数类型必填说明coord_states有序可迭代是广义坐标若给了speeds或全部状态right_hand_sideMatrix是F_1/F_2/F_3语义由其他参数决定speeds有序可迭代否广义速度提供则第一参数解释为坐标mass_matrixMatrix否M_2无coordinate_derivatives或M_3有coordinate_derivativesMatrix否提供即形式 [3]内容为Galg_con可迭代否代数约束行索引形式 [3] 下自动偏移为合并形式行号output_eqns字典否需跟踪的输出方程{名称: 表达式}coord_idxs/speed_idxs可迭代否状态中坐标/速度的索引bodies可迭代否Body/RigidBody对象集合loads可迭代否(point, force)或(frame, torque)元组集合主要属性与方法名称形态可用前提statesMatrix(o,1)总是可用coordinates/speedsMatrix初始化时能确定对应分量comb_explicit_rhsMatrix(o,1)形式 [1] 直接传入或先调用compute_explicit_form()comb_implicit_mat/comb_implicit_rhsMatrix(o,o)/Matrix(o,1)形式 [2] 传入或形式 [3] 自动拼接dyn_implicit_mat/dyn_implicit_rhsMatrix(m,m)/Matrix(m,1)仅形式 [3]kin_explicit_rhsMatrix(n,1)仅形式 [3]alg_conList初始化时提供返回合并形式行号bodies/loadsTuple初始化时提供dynamic_symbols()/constant_symbols()tuple总是可用compute_explicit_form()方法形式 [2] 或 [3]且未显式传入显式右端七、小结SymbolicSystem用一套紧凑的接口解决了“运动方程以什么格式交给数值代码”的问题三种形式覆盖了显式/隐式、运动学与动力学合并/分离的典型场景alg_con索引约定让 DAE 求解器能正确识别代数约束bodies/loads与output_eqns则为后续线性化、输出追踪预留了空间。其属性设计遵循“传入什么、自动推导什么、越界访问即报错”的明确边界配合 test_system.py 中三种形式等价性的完整断言三种形式构造的实例在comb_implicit_*与compute_explicit_form()之后输出一致可作为对接第三方数值求解器时行为契约的可靠依据。实际使用时建议从 3.6 节的最小构造方式起步确定方程属于哪种形式、核对alg_con的行号基准、区分“坐标速度”与“完整状态”两种传参风格即可获得一个可供数值分析代码消费的标准化系统描述。【免费下载链接】sympyA computer algebra system written in pure Python项目地址: https://gitcode.com/GitHub_Trending/sy/sympy创作声明:本文部分内容由AI辅助生成(AIGC),仅供参考
返回列表