物理/力学中的线性化¶

sympy.physics.mechanics 包括了在线性化生成的运动方程（EOM）关于一个工作点（也称为配平条件）的方法。请注意，这个工作点不一定是平衡位置，它只需要满足运动方程。

线性化是通过对操作点附近的运动方程进行一阶泰勒展开来实现的。当没有相关坐标或速度时，这仅仅是关于 \(q\) 和 \(u\) 的右侧的雅可比矩阵。然而，在存在约束的情况下，需要更加小心。这里提供的线性化方法正确处理了这些约束。

背景¶

在 sympy.physics.mechanics 中，我们假设所有系统都可以用以下一般形式表示：

\[\begin{split}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-m+k) \times 1}\end{split}\]

哪里

\[\begin{split}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{split}\]

在这种形式下，

\(f_{c}\) 表示配置约束方程
\(f_{v}\) 表示速度约束方程
\(f_{a}\) 表示加速度约束方程
\(f_{0}\) 和 \(f_{1}\) 构成了运动学微分方程
\(f_{2}\)、\(f_{3}\) 和 \(f_{4}\) 构成了动态微分方程
\(q\) 和 \(\dot{q}\) 是广义坐标及其导数
\(u\) 和 \(\dot{u}\) 是广义速度及其导数
\(r\) 是系统输入
\(\lambda\) 是拉格朗日乘数

这种广义形式保存在 Linearizer 类中，该类执行实际的线性化。KanesMethod 和 LagrangesMethod 对象都有使用 to_linearizer 类方法形成线性化器的方法。

一旦系统被强制转换为广义形式，就可以求解线性化的EOM。sympy.physics.mechanics 中提供的方法允许两种不同的线性化EOM形式：

\(M\), \(A\), 和 \(B\)

在这种形式中，强制矩阵被线性化为两个独立的矩阵 \(A\) 和 \(B\)。这是线性化运动方程的默认形式。得到的方程是：

\[ \begin{align}\begin{aligned}\begin{split}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}\end{split}\\\begin{split}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}\end{split}\end{aligned}\end{align} \]

哪里

\[\begin{split}M &\in \mathbb{R}^{(n+o+k) \times (n+o+k)}\\ A &\in \mathbb{R}^{(n+o+k) \times (n-l+o-m)}\\ B &\in \mathbb{R}^{(n+o+k) \times s}\end{split}\]

注意，\(q_i\) 和 \(u_i\) 只是独立的坐标和速度，而 \(q\) 和 \(u\) 包含独立和依赖的坐标和速度。

\(A\) 和 \(B\)

在这种形式中，线性化的运动方程被转换为显式的一阶形式，仅以独立坐标和速度表示。这种形式常用于稳定性分析或控制理论。得到的方程如下：

\[\begin{split}\begin{bmatrix} \delta \dot{q_i} \\ \delta \dot{u_i} \end{bmatrix} = A \begin{bmatrix} \delta q_i \\ \delta u_i \end{bmatrix} + B \begin{bmatrix} \delta r \end{bmatrix}\end{split}\]

哪里

\[ \begin{align}\begin{aligned}\begin{split}A &\in \mathbb{R}^{(n-l+o-m) \times (n-l+o-m)}\\ B &\in \mathbb{R}^{(n-l+o-m) \times s}\end{split}\\\begin{split}A &\in \mathbb{R}^{(n-l+o-m) \times (n-l+o-m)}\\ B &\in \mathbb{R}^{(n-l+o-m) \times s}\end{split}\end{aligned}\end{align} \]

要使用此表单，请在 linearize 类方法中设置 A_and_B=True。

线性化 Kane 方程¶

在初始化 KanesMethod 对象并使用 kanes_equations 类方法形成 \(F_r\) 和 \(F_r^*\) 之后，可以通过几种方式实现线性化。不同的方法将通过一个简单的摆系统来演示：

>>> from sympy import symbols, Matrix
>>> from sympy.physics.mechanics import *
>>> q1 = dynamicsymbols('q1')                     # Angle of pendulum
>>> u1 = dynamicsymbols('u1')                     # Angular velocity
>>> q1d = dynamicsymbols('q1', 1)
>>> L, m, t, g = symbols('L, m, t, g')

>>> # Compose world frame
>>> N = ReferenceFrame('N')
>>> pN = Point('N*')
>>> pN.set_vel(N, 0)

>>> # A.x is along the pendulum
>>> A = N.orientnew('A', 'axis', [q1, N.z])
>>> A.set_ang_vel(N, u1*N.z)

>>> # Locate point P relative to the origin N*
>>> P = pN.locatenew('P', L*A.x)
>>> vel_P = P.v2pt_theory(pN, N, A)
>>> pP = Particle('pP', P, m)

>>> # Create Kinematic Differential Equations
>>> kde = Matrix([q1d - u1])

>>> # Input the force resultant at P
>>> R = m*g*N.x

>>> # Solve for eom with kanes method
>>> KM = KanesMethod(N, q_ind=[q1], u_ind=[u1], kd_eqs=kde)
>>> fr, frstar = KM.kanes_equations([pP], [(P, R)])

1. 直接使用 `Linearizer` 类：¶

可以使用 to_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_B=True 来生成 \(A\) 和 \(B\) 形式：：

>>> A, B = linearizer.linearize(A_and_B=True)
>>> 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_B=True, op_point=op_point)
>>> A_op
Matrix([
[     0, 1],
[-g/L, 0]])

注意，通过将 msubs 应用于未使用 op_point 关键字参数生成的矩阵，可以获得相同的效果：

>>> assert msubs(A, op_point) == A_op

有时返回的矩阵可能不是最简化的形式。简化可以在事后进行，或者可以通过将 simplify kwarg 设置为 True 来使 Linearizer 对象在内部执行简化。

2. 使用 `linearize` 类方法：¶

KanesMethod 类的 linearize 方法提供了一个很好的包装器，它在内部调用 to_linearizer，执行线性化，并返回结果。请注意，上述 linearize 方法中可用的所有 kwargs 在这里也同样可用：

>>> A, B, inp_vec = KM.linearize(A_and_B=True, op_point=op_point, new_method=True)
>>> A
Matrix([
[     0, 1],
[-g/L, 0]])

额外的输出 inp_vec 是一个包含所有找到的 dynamicsymbols 的向量，这些符号不包含在广义坐标或速度向量中。这些被假定为系统的输入，形成上述背景中描述的 \(r\) 向量。在这个例子中没有输入，所以向量是空的：

>>> inp_vec
Matrix(0, 0, [])

线性化拉格朗日方程¶

拉格朗日方程的线性化过程与凯恩方程的过程非常相似。与之前一样，该过程将通过一个简单的摆系统来演示：

>>> # Redefine A and P in terms of q1d, not u1
>>> 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)

>>> # Solve for eom with Lagrange's method
>>> Lag = Lagrangian(N, pP)
>>> LM = LagrangesMethod(Lag, [q1], forcelist=[(P, R)], frame=N)
>>> lag_eqs = LM.form_lagranges_equations()

1. 直接使用 `Linearizer` 类：¶

可以通过 to_linearizer 类方法从 LagrangesMethod 对象形成一个 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_B=True, op_point=op_point)
>>> A
Matrix([
[     0, 1],
[-g/L, 0]])

2. 使用 `linearize` 类方法：¶

类似于 KanesMethod，LagrangesMethod 类也提供了一个 linearize 方法，作为一个很好的封装，它在内部调用 to_linearizer，执行线性化，并返回结果。与之前一样，唯一的区别是独立和从属坐标及速度也必须在调用中指定：

>>> A, B, inp_vec = LM.linearize(q_ind=[q1], qd_ind=[q1d], A_and_B=True, op_point=op_point)
>>> A
Matrix([
[     0, 1],
[-g/L, 0]])

潜在问题¶

虽然 Linearizer 类应该能够线性化所有系统，但可能会出现一些潜在问题。这些问题将在下面讨论，并提供一些解决它们的故障排除技巧。

1. 使用 `A_and_B=True` 进行符号线性化速度较慢¶

这可能是由于多种原因，但最可能的原因是，符号化地求解一个大型线性系统是一个昂贵的操作。指定一个工作点将减少表达式的大小并加快这一过程。然而，如果需要一个纯粹的符号化解决方案（例如，为了在以后的应用中处理许多工作点），一种解决方法是使用 A_and_B=False 进行评估，然后在应用工作点后手动求解：

>>> 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 和 M 中的符号越少，这个解决方案就会越快。因此，对于较大的表达式，推迟转换为 \(A\) 和 \(B\) 形式，直到大多数符号被替换为其数值，可能会对你有利。

2. 线性化形式具有 `nan`、`zoo` 或 `oo` 作为矩阵元素¶

这有两个可能的原因。第一个（也是你应该首先检查的）是某些依赖坐标的选项会在某些操作点导致奇点。以系统的方式进行坐标分区以避免这种情况超出了本指南的范围；更多信息请参见 [Blajer1994]。

造成这种情况的另一个潜在原因是，在进行操作点替换之前，矩阵可能没有处于最简化的形式。这种行为的一个简单例子是：

>>> 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 并设置 smart=True 将应用一种尝试避免这些条件的算法。然而，对于大型表达式，这会非常耗时。:

>>> msubs(expr, op_point, smart=True)
1

物理/力学中的线性化¶

背景¶

线性化 Kane 方程¶

1. 直接使用 Linearizer 类：¶

2. 使用 linearize 类方法：¶

线性化拉格朗日方程¶

1. 直接使用 Linearizer 类：¶

2. 使用 linearize 类方法：¶

潜在问题¶

1. 使用 A_and_B=True 进行符号线性化速度较慢¶

2. 线性化形式具有 nan、zoo 或 oo 作为矩阵元素¶

更多示例¶

1. 直接使用 `Linearizer` 类：¶

2. 使用 `linearize` 类方法：¶

1. 直接使用 `Linearizer` 类：¶

2. 使用 `linearize` 类方法：¶

1. 使用 `A_and_B=True` 进行符号线性化速度较慢¶

2. 线性化形式具有 `nan`、`zoo` 或 `oo` 作为矩阵元素¶