Autolev解析器#

介绍#

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

The parser has been built using the ANTLR framework and its main purpose is to help former users of Autolev to get familiarized with multibody dynamics in SymPy.

下面的部分将讨论解析器的细节，如用法、问题和未来的改进。有关Autolev和SymPy Mechanics的详细比较，您可能需要查看 SymPy Mechanics for Autolev Users guide .

使用#

我们首先从一个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()

注释的代码不是输出代码的一部分。print语句演示了如何获得与Autolev文件中的响应类似的响应。注意，我们需要使用SymPy函数，比如 .ang_vel_in() ， .dcm() 在很多情况下不像直接打印出变量 zero . 如果你对SymPy力学完全陌生 SymPy Mechanics for Autolev Users guide 导游应该会帮忙的。您可能还需要使用基本的SymPy简化和操作，例如 trigsimp() ， expand() ， evalf() 等获得类似于Autolev的输出。请参阅 SymPy Tutorial 想知道更多关于这些的事情。

高查斯#

  • 不要使用与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中的默认值。你也可以加上 np.deg2rad() 直接在密码里。

    对于输出代码（在分析 CODE 命令）作为解析器在 deg 单位在 INPUT 声明。

    这个 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 . 建议您在Symphy和Python中寻找替代方案。

  • 这个 REPRESENT 不支持命令。使用 MATRIX ， VECTOR 或 DYADIC 而是命令。Autolev4.1建议 REPRESENT 虽然仍然允许，但是解析器不解析它。

  • 不要使用该类型的变量声明 WO{{3}}RD{{2,4}} . 解析器只能处理一个变量名，后跟一对大括号和任意数量的 ' s、 如果您想实现以下目标，您必须手动声明所有案例 WO{{3}}RD{{2,4}} .

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

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

    这将计算 U1 ， U2 和 U3 在里面 E1 和 E2 . 最好使用这些命令的常规版本手动构造矩阵。

    %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 声明。解析器需要区分它们，以便将正确的参数传递给凯恩的方法对象。

    最好总是声明与坐标相对应的速度，并传入运动微分方程。解析器可以通过引入一些自己的伪变量来处理一些情况，但SymPy本身确实需要它们。

    还要注意，旧的Autolev声明像 VARIABLES U{{3}}' 也不支持。

    %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

  • Need to change me.dynamicsymbols._t to me.dynamicsymbols('t') for all occurrences of it in the Kane's equations. For example have a look at line 10 of this spring damper example. This equation is used in forming the Kane's equations so we need to change me.dynamicsymbols._t to me.dynamicsymbols('t') in this case.

    需要这样做的主要原因是PyDy要求显式地列出依赖时间的规范，而Autolev只需自己处理方程中的时间变量。

    问题是PyDy的系统类不接受 dynamicsymbols._t 作为指定的。参考问题 #396 . 这种变化实际上并不理想，因此将来应该找出更好的解决办法。

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

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

    另一方面，当一个声明的变量被分配给一个表达式时，解析器将针对该变量的表达式存储在字典中，这样就不会将它重新分配给一个完全不同的实体。这种限制是由于Python固有的特性以及它与Autolev等语言的不同之处。

    此外，Autolev似乎能够假设是使用变量还是在等式中指定变量的rhs表达式，即使没有显式表达式 RHS() 有时打电话。但是，要使解析器正常工作，最好使用 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] 不总是有效的，因此您可能需要查看用于求解返回值并正确索引的底层字典。

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

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

    2形状的惯性调整器 I_C_D>> = expr 但是，只有当C是一个实体时才设置惯量。如果C是粒子，那么 I_C_D>> = expr 简单地解析为 i_c_d = expr 和 i_c_d 就像一个规则变量。

    三。说到惰性吸气剂 (I_C_D>> 用于表达式或 INERTIA 命令），这些命令必须与 EXPRESS 命令将帧指定为SymPy需要此信息来计算惯性并矢。

    %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中，与标量名称一致的向量名将被覆盖。

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

    • 返回矩阵的一些方便形式的函数没有被解析。

    • 未分析设置。

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

    • 的已分析代码的字典索引 SOLVE 在许多情况下，命令是不恰当的。

    • 需要改变 dynamicsymbols._t 到 dynamicsymbols('t') 使PyDy模拟代码正常工作。

以下是一些其他的：

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

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

  • 相当于 TAYLOR 命令 .series() 不适用于 dynamicsymbols() .

  • 只有 DEPENDENT 当前正在分析约束。需要解析 AUXILIARY 约束也是如此。这应该尽快做，因为这不是很困难。

  • 现在没有能量和动量函数被解析。如果能让这些也能正常工作就好了。也许应该对SymPy做些改变。例如，SymPy没有等效于 NICHECK() .

  • 数值积分部分仅在 KANE 没有参数的命令。比如 KANE(F1,F2) 当前不工作。

  • 另外，PyDy数值模拟代码只适用于矩阵 ZERO = FR() + FRSTAR() 解决了。当矩阵中插入了其他一些方程时，它就不能很好地工作了。实现这一目标所面临的一个障碍是PyDy的系统类自动接受 forcing_full 和 mass_matrix_full 并在不给用户指定方程的灵活性的情况下解决这些问题。最好将此功能添加到系统类中。

未来改进#

1在线完成Dynamics#

The parser has been built by referring to and parsing codes from the Autolev Tutorial and the book Dynamics Online: Theory and Implementation Using Autolev. Basically, the process involved going through each of these codes, validating the parser results and improving the rules if required to make sure the codes parsed well.

这些代码的解析代码可以在GitLab上找到 here . 回购是私有的，因此需要请求访问。到现在为止，大多数代码到 动态在线 已被解析。

完成所有剩余代码， 2-10 ， 2-11 ， Ch4的其余部分 ， Ch5 和 Ch6 （不太重要）会使解析器更完整。

2解决问题#

要做的第二件事是着手解决上面在 Gotchas 和 Limitations and Issues 按优先级和易用性排序。其中许多都需要修改解析器代码，而有些则通过向SymPy添加一些功能来更好地修复。

三。切换到AST#

The parser is currently built using a kind of Concrete Syntax Tree (CST) using the ANTLR framework. It would be ideal to switch from a CST to an Abstract Syntax Tree (AST). This way, the parser code will be independent of the ANTLR grammar which makes it a lot more flexible. It would also be easier to make changes to the grammar and the rules of the parser.