Autolev 解析器

介绍

Autolev（现已被MotionGenesis取代）是一种用于符号多体动力学的领域特定语言。SymPy力学模块现在具有足够的功能和能力，可以成为一个功能齐全的符号动力学模块。此解析器通过利用SymPy的数学库和力学模块，将Autolev（4.1版本）代码解析为SymPy代码。

解析器是使用 ANTLR 框架构建的，其主要目的是帮助 Autolev 的以前用户熟悉 SymPy 中的多体动力学。

以下章节将讨论解析器的详细信息，如使用方法、注意事项、问题和未来改进。关于 Autolev 和 SymPy Mechanics 的详细比较，您可能需要查看 SymPy Mechanics for Autolev Users 指南。

用法

我们首先从一个 Autolev 代码文件开始。

让我们来看这个例子（注释 % 已包含在内，以显示 Autolev 的响应）：

% double_pendulum.al
%-------------------
MOTIONVARIABLES' Q{2}', U{2}'
CONSTANTS L,M,G
NEWTONIAN N
FRAMES A,B
SIMPROT(N, A, 3, Q1)
% -> N_A = [COS(Q1), -SIN(Q1), 0; SIN(Q1), COS(Q1), 0; 0, 0, 1]
SIMPROT(N, B, 3, Q2)
% -> N_B = [COS(Q2), -SIN(Q2), 0; SIN(Q2), COS(Q2), 0; 0, 0, 1]
W_A_N>=U1*N3>
% -> W_A_N> = U1*N3>
W_B_N>=U2*N3>
% -> W_B_N> = U2*N3>
POINT O
PARTICLES P,R
P_O_P> = L*A1>
% -> P_O_P> = L*A1>
P_P_R> = L*B1>
% -> P_P_R> = L*B1>
V_O_N> = 0>
% -> V_O_N> = 0>
V2PTS(N, A, O, P)
% -> V_P_N> = L*U1*A2>
V2PTS(N, B, P, R)
% -> V_R_N> = L*U1*A2> + L*U2*B2>
MASS P=M, R=M
Q1' = U1
Q2' = U2
GRAVITY(G*N1>)
% -> FORCE_P> = G*M*N1>
% -> FORCE_R> = G*M*N1>
ZERO = FR() + FRSTAR()
% -> ZERO[1] = -L*M*(2*G*SIN(Q1)+L*(U2^2*SIN(Q1-Q2)+2*U1'+COS(Q1-Q2)*U2'))
% -> ZERO[2] = -L*M*(G*SIN(Q2)-L*(U1^2*SIN(Q1-Q2)-U2'-COS(Q1-Q2)*U1'))
KANE()
INPUT M=1,G=9.81,L=1
INPUT Q1=.1,Q2=.2,U1=0,U2=0
INPUT TFINAL=10, INTEGSTP=.01
CODE DYNAMICS() some_filename.c

解析器可以如下使用:

>>> from sympy.parsing.autolev import parse_autolev
>>> sympy_code = parse_autolev(open('double_pendulum.al'), include_numeric=True)

# The include_pydy flag is False by default. Setting it to True will
# enable PyDy simulation code to be outputted if applicable.

>>> print(sympy_code)
import sympy.physics.mechanics as me
import sympy as sm
import math as m
import numpy as np

q1, q2, u1, u2 = me.dynamicsymbols('q1 q2 u1 u2')
q1d, q2d, u1d, u2d = me.dynamicsymbols('q1 q2 u1 u2', 1)
l, m, g=sm.symbols('l m g', real=True)
frame_n=me.ReferenceFrame('n')
frame_a=me.ReferenceFrame('a')
frame_b=me.ReferenceFrame('b')
frame_a.orient(frame_n, 'Axis', [q1, frame_n.z])
# print(frame_n.dcm(frame_a))
frame_b.orient(frame_n, 'Axis', [q2, frame_n.z])
# print(frame_n.dcm(frame_b))
frame_a.set_ang_vel(frame_n, u1*frame_n.z)
# print(frame_a.ang_vel_in(frame_n))
frame_b.set_ang_vel(frame_n, u2*frame_n.z)
# print(frame_b.ang_vel_in(frame_n))
point_o=me.Point('o')
particle_p=me.Particle('p', me.Point('p_pt'), sm.Symbol('m'))
particle_r=me.Particle('r', me.Point('r_pt'), sm.Symbol('m'))
particle_p.point.set_pos(point_o, l*frame_a.x)
# print(particle_p.point.pos_from(point_o))
particle_r.point.set_pos(particle_p.point, l*frame_b.x)
# print(particle_p.point.pos_from(particle_r.point))
point_o.set_vel(frame_n, 0)
# print(point_o.vel(frame_n))
particle_p.point.v2pt_theory(point_o,frame_n,frame_a)
# print(particle_p.point.vel(frame_n))
particle_r.point.v2pt_theory(particle_p.point,frame_n,frame_b)
# print(particle_r.point.vel(frame_n))
particle_p.mass = m
particle_r.mass = m
force_p = particle_p.mass*(g*frame_n.x)
# print(force_p)
force_r = particle_r.mass*(g*frame_n.x)
# print(force_r)
kd_eqs = [q1d - u1, q2d - u2]
forceList = [(particle_p.point,particle_p.mass*(g*frame_n.x)), (particle_r.point,particle_r.mass*(g*frame_n.x))]
kane = me.KanesMethod(frame_n, q_ind=[q1,q2], u_ind=[u1, u2], kd_eqs = kd_eqs)
fr, frstar = kane.kanes_equations([particle_p, particle_r], forceList)
zero = fr+frstar
# print(zero)
#---------PyDy code for integration----------
from pydy.system import System
sys = System(kane, constants = {l:1, m:1, g:9.81},
specifieds={},
initial_conditions={q1:.1, q2:.2, u1:0, u2:0},
times = np.linspace(0.0, 10, 10/.01))

y=sys.integrate()

注释的代码不是输出代码的一部分。打印语句展示了如何获得类似于Autolev文件中的响应。请注意，在许多情况下我们需要使用SymPy函数，如``.ang_vel_in()``、.dcm()``等，而不是像``zero``那样直接打印变量。如果你对SymPy力学完全陌生，:ref:`SymPy Mechanics for Autolev Users指南 <sympy_mechanics_for_autolev_users>` 应该会有所帮助。你可能还需要使用基本的SymPy简化与操作，如``trigsimp()、expand()、``evalf()``等，以获得类似于Autolev的输出。请参阅 SymPy教程 以了解更多信息。

陷阱

  • 不要使用与Python保留字冲突的变量名。这是一个违反此规则的例子：

    %Autolev Code
%------------
LAMBDA = EIG(M)

    #SymPy Code
#----------
lambda = sm.Matrix([i.evalf() for i in (m).eigenvals().keys()])

  • 确保向量和标量的名称不同。Autolev 对这些处理方式不同，但在 Python 中这些会被覆盖。解析器目前允许刚体和标量/向量的名称相同，但不允许标量和向量之间的名称相同。这可能在将来需要改变。

    %Autolev Code
%------------
VARIABLES X,Y
FRAMES A
A> = X*A1> + Y*A2>
A = X+Y

    #SymPy Code
#----------
x, y = me.dynamicsymbols('x y')
frame_a = me.ReferenceFrame('a')
a = x*frame_a.x + y*frame_a.y
a = x + y
# Note how frame_a is named differently so it doesn't cause a problem.
# On the other hand, 'a' gets rewritten from a scalar to a vector.
# This should be changed in the future.

  • 当处理由函数返回的矩阵时，必须检查值的顺序，因为它们可能与在Autolev中的顺序不同。特别是对于特征值和特征向量，这种情况尤为明显。

    %Autolev Code
%------------
EIG(M, E1, E2)
% -> [5; 14; 13]
E2ROW = ROWS(E2, 1)
EIGVEC> = VECTOR(A, E2ROW)

    #SymPy Code
#----------
e1 = sm.Matrix([i.evalf() for i in m.eigenvals().keys()])
# sm.Matrix([5;13;14]) different order
e2 = sm.Matrix([i[2][0].evalf() for i in m.eigenvects()]).reshape(m.shape[0], m.shape[1])
e2row = e2.row(0)
# This result depends on the order of the vectors in the eigenvecs.
eigenvec = e2row[0]*a.x + e2row[1]*a.y + e2row[2]*a.y

  • 在使用 EVALUATE 时，使用类似 90*UNITS(deg,rad) 的角度替换，因为在 SymPy 中弧度是默认的。你也可以直接在 SymPy 代码中添加 np.deg2rad()。

    对于输出代码（在解析 CODE 命令时生成），不需要这样做，因为当 INPUT 声明中给出 deg 单位时，解析器会处理这个问题。

    另一方面，DEGREES 设置只在某些情况下有效，例如在 SIMPROT 中，当需要一个角度时。

    %Autolev Code
%------------
A> = Q1*A1> + Q2*A2>
B> = EVALUATE(A>, Q1:30*UNITS(DEG,RAD))

    #SymPy Code
#----------
a = q1*a.frame_a.x + q2*frame_a.y
b = a.subs({q1:30*0.0174533})
# b = a.subs({q1:np.deg2rad(30)}

  • 大多数 Autolev 设置尚未被解析，对解析器没有影响。唯一有点作用的是 COMPLEX 和 DEGREES。建议在 SymPy 和 Python 中寻找这些的替代方案。

  • REPRESENT 命令不受支持。请改用 MATRIX、VECTOR 或 DYADIC 命令。Autolev 4.1 也建议使用这些命令而不是 REPRESENT，尽管仍然允许使用，但解析器不会解析它。

  • 不要使用 WO{3}RD{2,4} 类型的变量声明。解析器只能处理一个变量名后面跟着一对花括号和任意数量的 ' 。如果你想实现类似 WO{3}RD{2,4} 的效果，你必须手动声明所有的情况。

  • 解析器可以处理大多数命令的常规版本，但在大多数情况下可能无法正确解析带有矩阵参数的函数。例如：

    M=COEF([E1;E2],[U1,U2,U3])

    这将计算 E1 和 E2 中 U1、U2 和 U3 的系数。使用这些命令的常规版本手动构建矩阵是更可取的。

    %Autolev Code
%------------
% COEF([E1;E2],[U1,U2,U3])
M = [COEF(E1,U1),COEF(E1,U2),COEF(E1,U3) &
    ;COEF(E2,U1),COEF(E2,U2),COEF(E2,U3)]

  • MOTIONVARIABLE 声明必须用于广义坐标和速度，以及所有其他变量必须在常规的 VARIABLE 声明中声明。解析器需要这一点来区分它们，以便将正确的参数传递给 Kane 方法对象。

    同样建议始终声明与坐标对应的速度，并传入运动学微分方程。解析器能够在某些情况下处理未声明的情况，方法是引入一些自己的虚拟变量，但SymPy本身确实需要这些声明。

    另请注意，像 VARIABLES U{3}' 这样的旧版 Autolev 声明也不受支持。

    %Autolev Code
%------------
MOTIONVARIABLES' Q{2}', U{2}'
% ----- OTHER LINES ----
Q1' = U1
Q2' = U2
----- OTHER LINES ----
ZERO = FR() + FRSTAR()

    #SymPy Code
#----------
q1, q2, u1, u2 = me.dynamicsymbols('q1 q2 u1 u2')
q1d, q2d, u1d, u2d = me.dynamicsymbols('q1 q2 u1 u2', 1)

# ------- other lines -------

kd_eqs = [q1d - u1, q2d - u2]
kane = me.KanesMethod(frame_n, q_ind=[q1,q2], u_ind=[u1, u2], kd_eqs = kd_eqs)
fr, frstar = kane.kanes_equations([particle_p, particle_r], forceList)
zero = fr+frstar

  • 需要将 Kane 方程中所有出现的 me.dynamicsymbols._t 更改为 me.dynamicsymbols('t')。例如，查看此 弹簧阻尼器示例 的第 10 行。此方程用于形成 Kane 方程，因此在这种情况下我们需要将 me.dynamicsymbols._t 更改为 me.dynamicsymbols('t')。

    需要这样做的一个主要原因是，PyDy 需要明确地列出时间相关的指定项，而 Autolev 则自行处理方程中的游离时间变量。

    问题是 PyDy 的 System 类不接受 dynamicsymbols._t 作为指定的参数。请参阅问题 #396。这个改变实际上并不理想，所以未来应该找到一个更好的解决方案。

  • 解析器通过解析 Autolev 代码中的变量声明来创建 SymPy symbols 和 dynamicsymbols。

    对于直接初始化的中间表达式，解析器不会创建 SymPy 符号。它只是将它们分配给表达式。

    另一方面，当一个声明的变量被赋值给一个表达式时，解析器会将该表达式存储在字典中以防止变量被重新赋值给一个完全不同的实体。这种约束是由于Python的固有特性及其与Autolev等语言的区别所导致的。

    此外，Autolev 似乎能够在某些情况下假设是否在方程中使用变量或该变量已被赋值的右侧表达式，即使在没有显式 RHS() 调用的情况下也是如此。然而，为了使解析器正确工作，最好在任何需要使用变量右侧表达式的地方使用 RHS()。

    %Autolev Code
%------------
VARIABLES X, Y
E = X + Y
X = 2*Y

RHS_X = RHS(X)

I1 = X
I2 = Y
I3 = X + Y

INERTIA B,I1,I2,I3
% -> I_B_BO>> = I1*B1>*B1> + I2*B2>*B2> + I3*B3>*B3>

    #SymPy Code
#----------
x,y = me.dynamicsymbols('x y')
e = x + y  # No symbol is made out of 'e'

# an entry like {x:2*y} is stored in an rhs dictionary

rhs_x = 2*y

i1 = x  # again these are not made into SymPy symbols
i2 = y
i3 = x + y

body_b.inertia = (me.inertia(body_b_f, i1, i2, i3), b_cm)
# This prints as:
# x*b_f.x*b_f.x + y*b_f.y*b_f.y + (x+y)*b_f.z*b_f.z
# while Autolev's output has I1,I2 and I3 in it.
# Autolev however seems to know when to use the RHS of I1,I2 and I3
# based on the context.

  • 这是 SOLVE 命令的解析方式：

    %Autolev Code
%------------
SOLVE(ZERO,X,Y)
A = RHS(X)*2 + RHS(Y)

    #SymPy Code
#----------
print(sm.solve(zero,x,y))
# Behind the scenes the rhs of x
# is set to sm.solve(zero,x,y)[x].
a = sm.solve(zero,x,y)[x]*2 + sm.solve(zero,x,y)[y]

    像 [x] 和 [y] 这样的索引并不总是有效，因此你可能需要查看 solve 返回的基础字典，并正确地对其进行索引。

  • 在解析器的上下文中，惯性声明和惯性函数的工作方式有些不同。这可能一开始很难理解，但由于SymPy和Autolev之间的差异，必须这样做以弥合差距。以下是关于它们的一些要点：

    惯性声明 (INERTIA B,I1,I2,I3) 设置刚体的惯性。

    2. Inertia setters of the form I_C_D>> = expr however, set the inertias only when C is a body. If C is a particle then I_C_D>> = expr simply parses to i_c_d = expr and i_c_d acts like a regular variable.

    3. When it comes to inertia getters (I_C_D>> used in an expression or INERTIA commands), these MUST be used with the EXPRESS command to specify the frame as SymPy needs this information to compute the inertia dyadic.

    %Autolev Code
%------------
INERTIA B,I1,I2,I3
I_B_BO>> = X*A1>*A1> + Y*A2>*A2>  % Parser will set the inertia of B
I_P_Q>> = X*A1>*A1> + Y^2*A2>*A2> % Parser just parses it as i_p_q = expr

E1 = 2*EXPRESS(I_B_O>>,A)
E2 =  I_P_Q>>
E3 = EXPRESS(I_P_O>>,A)
E4 = EXPRESS(INERTIA(O),A)

% In E1 we are using the EXPRESS command with I_B_O>> which makes
% the parser and SymPy compute the inertia of Body B about point O.

% In E2 we are just using the dyadic object I_P_Q>> (as I_P_Q>> = expr
% doesn't act as a setter) defined above and not asking the parser
% or SymPy to compute anything.

% E3 asks the parser to compute the inertia of P about point O.

% E4 asks the parser to compute the inertias of all bodies wrt about O.

  • 在物体的惯性声明中，如果惯性是围绕质心以外的点设置的，则需要确保该点的位置向量设置器和质心出现在惯性声明之前，否则 SymPy 将抛出错误。

    %Autolev Code
%------------
P_SO_O> = X*A1>
INERTIA S_(O) I1,I2,I3

  • 请注意，并非所有 Autolev 命令都已实现。解析器现在涵盖了它们基本形式中的重要命令。如果您不确定某个命令是否包含在内，请查看源代码中的 this file。搜索“<command>”以验证这一点。查看特定命令的代码也会让您了解它预期的工作形式。

限制与问题

  • 很多问题已经在 Gotchas 部分讨论过了。其中一些包括：

    • 与标量名称相同的向量名称在Python中会被覆盖。

    • 一些方便的变量声明未被解析。

    • 一些返回矩阵的函数的便捷形式无法解析。

    • 设置未被解析。

    • 符号和右侧表达式在Python中的工作方式非常不同，这可能会导致不理想的结果。

    • 在许多情况下，SOLVE 命令的解析代码的字典索引不正确。

    • 需要将 dynamicsymbols._t 更改为 dynamicsymbols('t') 以使 PyDy 模拟代码正常工作。

以下是一些其他的：

  • 特征向量似乎没有按预期工作。在许多情况下，Autolev 和 SymPy 中的值并不相同。

  • 块矩阵不会被解析器解析。实际上，在 SymPy 中进行更改以允许矩阵接受其他矩阵作为参数会更容易。

  • SymPy 中 TAYLOR 命令的等价物 .series() 不适用于 dynamicsymbols()。

  • 目前仅解析 DEPENDENT 约束。还需要解析 AUXILIARY 约束。这应该很快完成，因为它并不非常困难。

  • 目前没有任何能量和动量函数被解析。如果能把这些功能也实现就更好了。可能需要对SymPy做一些改动。例如，SymPy没有与 NICHECK() 等效的函数。

  • 数值积分部分仅在 KANE 命令无参数的情况下正常工作。类似 KANE(F1,F2) 的功能目前无法使用。

  • 此外，PyDy 数值模拟代码仅适用于求解矩阵 ZERO = FR() + FRSTAR() 的情况。当矩阵中插入其他方程时，它不能很好地工作。实现这一目标面临的一个障碍是，PyDy 的 System 类自动接收 forcing_full 和 mass_matrix_full 并求解它们，而没有给用户指定方程的灵活性。为 System 类添加此功能将会很有帮助。

未来改进

1. 完成在线动态

解析器是通过参考和解析 Autolev 教程 以及书籍 Dynamics Online: Theory and Implementation Using Autolev 中的代码构建的。基本上，这个过程涉及逐一检查这些代码，验证解析器的结果，并在必要时改进规则，以确保代码解析良好。

这些解析后的代码可以在 GitLab 上找到 这里 。该仓库是私有的，因此需要请求访问权限。截至目前，大部分代码已解析至《Dynamics Online》的第4章。

完成书中剩余的所有代码（即 2-10、2-11、第4章的其余部分、第5章 和 *第6章*（不太重要））将使解析器更加完整。

2. 修复问题

第二件事是要按照优先级和难易程度来解决 陷阱 和 限制与问题 部分中描述的问题。其中许多问题需要修改解析器代码，而有些问题则通过向 SymPy 添加一些功能来更好地解决。

3. 切换到AST

解析器目前使用的是基于 ANTLR 框架的某种具体语法树（CST）构建的。理想情况下，应该从 CST 切换到抽象语法树（AST）。这样，解析器代码将独立于 ANTLR 语法，从而使其更加灵活。修改语法和解析器规则也会变得更加容易。