本
文
摘
要
单摆
2.1 介绍
本章的目标是了解单摆的动力学。为什么我们要选择单摆来作为我们的研究对象呢？因为在某种程度上来说是多连杆机器人动力学是多个单摆动力学的耦合。此外，单摆的动力学分析内容丰富，涵盖了本书涉及的非线性动力学中的大部分基础概念，也恰好能够在很短的篇幅内将其讲清楚。
图2.1 单摆如附录所述，我们可以利用拉格朗日方程来推导单摆的运动方程：
ml2θ¨(t)+mglsinθ(t)=Qml^2\ddot{\theta}(t) + mgl \sin\theta(t)=Q
等式的右边是广义力 QQ ，包括摩擦产生的阻尼转矩和输入的控制力矩 u(t)u(t) ，其表达式为：
Q=−bθ˙(t)+u(t)Q=-b\dot{\theta}(t)+u(t)
2.2 恒定的扭矩下的非线性动力学
让我们首先考虑恒定扭矩驱动单摆的动力学方程：
ml2θ¨(t)+bθ˙(t)+mglsinθ(t)=u0ml^2\ddot{\theta}(t) + b\dot{\theta}(t)+ mgl \sin\theta(t)=u_0
实例2.1 Python的单摆 你可以尝试在Drake使用单摆实例
上述运动方程是相对简单的微分方程。但如果给定系统的初始值 θ(0){\theta}(0) 和 θ˙(0)\dot{\theta}(0) ，是否能够通过积分获得 θ(t){\theta}(t) 呢？虽然大多时候答案是肯定的，但是即使在比如说 b=u=0b=u=0 这样的简单的情况下，积分过程也会涉及第一类椭圆积分，而无法通过直觉和简单计算直接获得结果。 这与线性系统形成了鲜明的对比，线性系统能够显式地直接积分。举一个简单的线性系统的例子，我们有：
q˙=aq→q(t)=q(0)eat\dot{q}=aq \rightarrow q(t)=q(0)e^{at}
从积分结果我们可以马上知道系统的长期行为：
如果 2.2.1 过阻尼摆让我们先从学习一个特例开始。直观上来说当 bθ˙≫ml2θ¨b\dot{\theta} \gg ml^2\ddot{\theta} 时，也就是说 bgl≫ml2b\sqrt{\frac{g}{l}} \gg ml^2 时（通过单位分析得到 gl\sqrt{\frac{g}{l}} 是系统的本征频率），过阻尼情况就出现了，一个典型的实例就是单摆在糖浆中运动，此时阻尼项 bθ˙b\dot{\theta} 在加速度组成中占主导地位，那么我们的动力学方程可近似成：
ml2θ¨+bθ˙≈bθ˙=u0−mglsinθml^2\ddot{\theta}+b\dot{\theta} \approx b\dot{\theta}=u_0-mgl \sin\theta
换句话讲，在强过阻尼情况下，系统将近似变成一阶的，这是过阻尼系统的通用特征，常发生在极低雷诺数的流体中。 在这里我们忽略 θ\theta 的周期重复性，我们可以将系统的运动方程改写成
bx˙=u0−mglsinxb\dot{x}=u_0-mgl \sin x
我们的目标是了解系统的长期行为：即给定初始值 x(0)x(0) 找到 x(∞)x(\infty) 。让我们在 u0=0u_0=0 的情况下开始绘制 x˙\dot{x} 对 xx 的关系图：
图2.2 相图首先要注意的是，该系统具有多个_不动点_或者说_稳定状态_，即是对应于 x˙=0\dot{x}=0 时的状态。在这个简单的例子中，不动点 *** 为 x∗={...,−π,0,π,2π,...}x^*=\{...,-\pi,0,\pi,2\pi,...\} 。当系统处于任一不动点状态时，它永远不会离开不动点状态。也就是说，如果系统的初始状态是位于某个不动点，那么 x(∞)x(\infty) 也将会保持在同一不动点。 接下来让我们来探讨系统在不动点附近的行为。考虑上图中的不动点 x∗=πx^*=\pi ，如果系统在不动点的右边启动，那么此时 x˙\dot{x} 是正的，那么系统将远离不动点 x∗=πx^*=\pi ；如果系统在不动点左边启动，那么此时 x˙\dot{x} 是负的，那么系统将向负方向移动，远离不动点 x∗=πx^*=\pi 。我们称类似于 x∗=πx^*=\pi 的不动点是_不稳定的不动点_。如果我们再关注一下不动点 x∗=0x^*=0 ，那么情况是截然相反的：轨迹从不动点附近出发向右或向左移动最终回到不动点。我们把 x∗=0x^*=0 称为_局部稳定的_不动点。更具体地说，我们稳定性进行分类（其中 ϵ\epsilon 用于表示任意_小的_标量）：
李雅普诺夫意义上的局部稳定：对于任意的 0">ϵ>0\epsilon>0 的球内出发，即 例如，我们可以看到任意初始条件 x(0)∈(−π,π)x(0) \in (-\pi,\pi) 将导致 limt→∞x(t)=0\lim_{t \to \infty}x(t) = 0 ，这个区域 (−π,π)(-\pi,\pi) 被称为不动点 x∗=0x^*=0 的吸引域。两个不动点的吸引域不能重叠，分离两个吸引域的流形被称为_分界线_。不稳定的不动点 x∗={...,−π,π,3π,...}x^*=\{...,-\pi,\pi,3\pi,...\} 形成了由稳定不动点形成的吸引域之间的_分界线_。 正如图像显示，一阶一维系统的行为是相对受限的：单调逼近不动点或单调走向 ±∞\pm \infty ，而没有其他的可能性，振荡是不可能出现的。对于一阶非线性系统 （不仅仅是抛物线）来说，图形分析是一种出色的分析工具，正如下面的例子所示：实例2.2 非线性autapse 考虑以下系统：
x˙+x=tanh(wx)\dot{x}+x=\tanh(wx)
下图绘制两个不同 ww 值所对应的系统相图，注意到对于小 zz ，有 tanh(z)≈z\tanh(z)\approx z 。当 w≤1w \le 1 ，该系统仅具有单个不动点。当 1">w>1w > 1 ，该系统有三个不动点： 两个稳定和一个不稳定。
图2.4 非线性autapse这些方程并不是我们任意构造的，实际上其中一个是最简单的神经网络模型，另一个是最简单的持久性存储器模型[2]。在方程中， xx 代表单个神经元的放电模型，它有一个连接到本身的反馈。 tanh\tanh 是神经元的激活函数， ww 是突触反馈的权重。 可以自己尝试代码。另外，您还可以尝试寻找一种 LSTM单元的方程，看看你是否能找出答案!
介绍最后一条术语。在神经元和许多动力系统的例子中，可以通过单一参数对动力学参数化。当我们改变 ww ，系统的不动点会发生变化。事实上，如果我们增加 ww 到 w=1w =1 ，戏剧性的一幕发生了，系统的不动点的数量从一个变成三个，这就是所谓的_分岔_。这个特定的分支被称为pitchfork分岔。我们经常以参数为横坐标，系统的不动点为纵坐标来表示分岔图，用实线表示稳定不动点，虚线表示不稳定不动点，如图所示：
图2.5 非线性autapse的分岔图当改变单摆系统的恒定扭矩输入 u0u_0 时，我们也能得到一个分岔图。 最后，让我们回到原来以 θ\theta 为变量的方程，而不是 xx 。我们需要唯一注意的一点是：因为角度的循环性，这个系统看起来像是振荡的。但事实上，图形分析的结果告诉我们当 mgl">|u0|>mgl|u_0|>mgl 时的不稳定性。 递归神经网络是另一个很好的例子，在后面的Hopfield网络练习中，你也能发现不动点、吸引域和分叉这些概念的应用。
2.2.2 零扭矩无阻尼摆
再考虑系统：
ml2θ¨=u0−mglsinθ−bθ˙ml^2\ddot{\theta} = u_0 - mgl \sin \theta - b\dot{\theta}
这一次我们使 b=0b=0 ，系统的动力学就是二阶的了。我们总能把任何二阶系统分解成耦合的一阶系统，考虑不依赖于时间的一般二阶系统，其方程为：
q¨=f(q,q˙,u)\ddot{q}=f(q,\dot{q},u)
该二阶系统等同于二个一阶系统的耦合：
x˙1=x2x˙2=f(x1,x2,u)\begin{align} \dot{x}_1&=x_2 \\ \dot{x}_2&=f(x_1, x_2, u) \end{align}
其中 x1=qx_1= q 和 x2=q˙x_2=\dot{q} 。这个系统图形化的表示就不是一条线了，而是一种由矢量场的矢量 [x˙1T,x˙2T][\dot{x}_1^T,\dot{x}_2^T] 绘制的域 (x1,x2)(x_1,x_2) ，此向量场被称为系统的_相图_ 。 在本节中，我们只讨论 u0=0u_0 = 0 的简单情况，让我们先来绘制相图，首先沿 θ\theta 轴，矢量场的 xx 分量在这里为零，矢量场的 yy 分量为 −glsinθ-\frac{g}{l}\sin \theta 。正如预期的那样，不动点为 ±π...\pm \pi ... ，现在画出向量场的其余部分，你能告诉我哪个不动点是稳定吗？其中一些是李雅普诺夫意义稳定，但没有一个是渐近稳定的。
图2.6 二维相图轨道计算
你可能想知道我们是如何在上图中绘制黑色轮廓线的：我们可以通过数值模拟系统来获得它们，也可以很容易地以封闭形式获得。直接对运动方程积分是困难的，但至少在 u0=0u_{0}=0 时，一些额外的物理世界的见解可以被我们利用。单摆的动能 TT 和势能 UU 如下：
T=12Iθ˙2,U=−mglcos(θ)T=\frac{1}{2} I \dot{\theta}^{2}, \quad U=-m g l \cos (\theta)
其中 I=ml2I=m l^{2} 。系统的总能量是 E(θ,θ˙)=T(θ˙)+U(θ)E(\theta, \dot{\theta})=T(\dot{\theta})+U(\theta) 。无阻尼单摆是一个能量守恒系统，即系统总能量在轨迹上保持恒定。利用能量守恒，我们有：
E(θ(t),θ˙(t))=E(θ(0),θ˙(0))=E012Iθ˙2(t)−mglcos(θ(t))=E0θ˙(t)=±2I[E0+mglcos(θ(t))]\begin{array}{l} E(\theta(t), \dot{\theta}(t))=E(\theta(0), \dot{\theta}(0))=E_{0} \\ \quad \frac{1}{2} I \dot{\theta}^{2}(t)-m g l \cos (\theta(t))=E_{0} \\ \dot{\theta}(t)=\pm \sqrt{\frac{2}{I}\left[E_{0}+m g l \cos (\theta(t))\right]} \end{array}
利用上述方程，如果我们知道 θ\theta ，就可以确定 θ˙\dot{\theta} 两个可能的解，也就能表达图中黑色轮廓线的所有特性。当 \cos \left(\theta_{\max }\right)">cos(θ)>cos(θmax)\cos (\theta)>\cos \left(\theta_{\max }\right)
虽然控制器看起来有点奇怪，但它被选中的原因也很简单，因为事实证明由此产生的 “误差动力学”的形式也很简洁：
E~˙=−kθ˙2E~\dot{\tilde{E}}= -k\dot{\theta}^2\tilde{E}
对系统进行图形分析，如果绘制任意固定 θ˙\dot{\theta} 对应的 E~˙\dot{\tilde{E}} 与 E~{\tilde{E}} 的关系图，你会发现，它是一条分开水平轴的直线，也就将意味着系统满足指数收敛： E~→0\tilde{E} \to 0 ，这适用于除外 θ˙=0\dot{\theta} = 0 任何 θ˙\dot{\theta} 。因为 θ˙=0\dot{\theta} = 0 时单摆不会从下面的不动点往上摆动，但如果你轻推一下系统，单摆便会开始蓄积能量，便可以往上摆动了。最重要的属性是当 E^d">E>EdE > E^d ，我们应该像负阻尼一样增加能量。即使控制输入的扭矩是有界的，收敛也很容易被保持。 这是一个将系统的所有轨迹推向不稳定的不动点的非线性控制器，但它能使不稳定平衡点变成局部稳定吗？如果只有这个控制器，不动点是_吸引的，但不稳定_ ，就像我们上面的例子所示。一个小扰动可能会导致系统为了再次
实例 2.4 能量改造的摆
不妨会用一分钟的时间去尝试用不同能量改造控制器的摆
请确保花点时间看看这些实例代码，注意那些用来切换到平衡控制器的随机阈值。我们将在李雅普诺夫方法那一章介绍令人满意的求解方法。
这个控制器有几点非常好。它不仅很简单，而且在实际应用中对模型参数 m ， l 和 g 的估计误差非常地稳健。事实上，模型参数唯一影响控制方程的地方是在计算 E=\frac{1}{2}ml^2\dot{\theta}^2-mgl(1+ \cos \theta) 时，影响稳定性只在 \tilde{E} = 0 时对应的轨道的位置 \frac{1}{2}l \dot{\theta}^2 = g(1+ \cos \theta) 时，而整个过程不依赖于估计的质量 m ，只线性地依赖于估计的长度 l 和重力常数 g （如果不能准确测量长度单摆，那么更换更好的测量装置显然是一笔划得来的投资）。我们也可以开发多种优化算法来实现单摆的摆动与平衡，但找到一个对模型参数不敏感的模型需要大量的工作，如果原系统中有阻尼，我们也可以使用 u = -k\dot{\theta} \tilde{E} + b \dot{\theta} 来抵消阻尼。和前面一样，控制器对有误差的阻尼比估计是很稳健的。
2.4 练习
练习2.4.1 （图形分析）
图2.8 一阶系统的多重均衡考虑一阶系统：
\dot{x}=\left\{\begin{array}{lll} -x^{5}+2 x^{3}-x & \text { if } & x \leq 1 \\ 0 & \text { if } & 1<x \leq 2 \\ -x+2 & \text { if } & x>2 \end{array}\right.
其动力学如上图所示。请注意，包括 x^{*}=-1 和 x^{*}=0 在内的闭合区间 [1,2] 内的所有点都是这个系统的不动点。对于不动点 *** ，请你确定哪些是不稳定的、李雅普诺夫稳定的、渐近稳定的还是指数稳定的。
练习2.4.2 吸引域和分岔图
考虑具有多项式动力学的一阶系统：
\dot{x}=f(x)=-x^{3}+4 x^{2}+11 x-30
a.绘制函数图 f(x) ，并确定所有平衡点。（请随意选用工具来绘制此图以作验证）
b.对于每个不动点，确定它是李雅普诺夫稳定的还是不稳定的。对于每一个稳定不动点，确定吸引域。
c.考虑系统动力学中的附加项 w ： \dot{x}=f(x)+w 。随着 w 从 0 增加到 \infty ，确定系统有多少个稳定和不稳定的平衡点。绘制非负值 w 的分叉图来支持你的答案。
练习2.4.3 枚举不稳定平衡点
给定系统动力学 \dot{x}=f(x) ，具有标量状态 x 。动力学 f(x) 未知，但我们得到了两条信息：
函数 f(x) 是连续的;该系统恰好有三个李雅普诺夫稳定的平衡点。确定系统可能具有的不稳定不动点的最小值 n_{\min } 和最大值 n_{\max } 。绘制与 n_{\min } 和 n_{\max } 对应的 f_{\min }(x) 和 f_{\max }(x) 两个函数的草图来支持你的答案验证上述要求。
练习2.4.4 吸引性与稳定性
图2.9 二维系统相图上图显示了下面公式对应的动力学系统的相图：
\begin{array}{l} \dot{x}_{1}=x_{1}(1-|\mathbf{x}|)+x_{2} \frac{x_{1}-|\mathbf{x}|}{2|\mathbf{x}|} \\ \dot{x}_{2}=x_{2}(1-|\mathbf{x}|)-x_{1} \frac{x_{1}-|\mathbf{x}|}{2|\mathbf{x}|} \end{array}
其中 |\mathbf{x}|=\sqrt{x_{1}^{2}+x_{2}^{2}} 。为了帮助您分析这个系统，我们提供了这个python笔记本。花点时间理解其中的代码，然后回答以下问题。
a.找到该系统的所有平衡点 （无需查看上面描述的状态空间部分之外的部分）。使用笔记本仔细检查你的答案： 确定系统的不动点并将不动点设置为初始条件，模拟系统的演变。
b.确定你在前一点确定的不动点是否有吸引力或李雅普诺夫稳定的。简要解释你的推理，你也可以在书面回答中包含一些用笔记本生成的图。
c.这个动力系统与我们在本章中分析的某一个系统密切相关。你能猜出是哪一个吗？两者之间有什么联系？思考题： 用数学推导来支持你的主张。
练习2.4.5 振动基座单摆
考虑一个受驱动的单摆，它的底座（杆的旋转中心）根据谐波定律做水平振荡 h \sin (\omega t) ，其中 h 和 \omega 分别表示振荡的幅度和振荡的频率，其运动方程是：
m l^{2} \ddot{\theta}+m g l \sin \theta=m l \omega^{2} h \sin (\omega t) \cos \theta+u
（如果这个方程对你来说很模糊，试着用拉格朗日方法来推导它； 但是要小心，这个系统的动能 T(\theta, \dot{\theta}, t) 显式地取决于时间）我们的目标是设计一个时间相关的控制律 u=\pi(\theta, \theta, t) ，使得单摆以恒定的 \dot{\theta}=1 速度旋转。
a.求出使不动点唯一稳定且位于 \dot{\theta}=1 所需的闭环动力学 \ddot{\theta}=f(\dot{\theta}) 。
b.使用反馈抵消来设计控制律 u=\pi(\theta, \dot{\theta}, t) ，使得闭环动力学与在前一点选择的一致。
c.在这个Python笔记本里，我们建立了一个模拟环境来测试您的控制律。试着理解代码的结构： 这个工作流程相当普遍，它可能是你未来使用Drake完成项目的一个很好的起点。在专用单元中实施您的控制律，并使用笔记本末尾的动画和绘图检查您的工作。
练习2.4.6 Hopfield网络
考虑以下动力系统：
\dot{x}=A^{T} \operatorname{softmax}(\beta A x)-x
其中 A \in \mathbb{R}^{m \times n} 和 \beta \in \mathbb{R} ，分别是常数矩阵和标量，softmax函数由 \operatorname{softmax}(z)_{i}=\frac{e^{z_{i}}}{\sum_{k} e^{z_{k}}} 给出。
a.当 \beta=5 和 A=\left[\begin{array}{cc}1 & 0 \\ 0 & 1 \\ -1 & -1\end{array}\right] 时，绘制系统的相图。
b.通过观测相图来确定上述系统的不动点，稳定和不稳定的不动点是哪些？
c.矩阵 A 对应的稳定不动点是？
d.在Python笔记本中，使用Hopfield网络实现图像恢复算法。修改上述动力系统以记忆MNIST数据集的一个子集。使用损坏的图像作为初始条件，模拟动力学系统恢复图像的过程。
练习2.4.7 神经ODE图形分析
我们最近看到了一种使用神经网络来描述动力学函数的趋势。使用网络来表示微分方程可能非常复杂，也很难进行严格的分析。然而在一维系统这种特殊情况下，我们可以使用图形分析来分析我们的系统。在本练习中，您将尝试从数据中学习阻尼摆的动力学。本练习有一个编码部分，您将在其中学习在pytorch中训练神经网络的基础知识。你也可以利用已有知识，对代表动力学函数的训练后的网络进行图形分析。
a.在这个python笔记本中，从数据中学习阻尼单摆的动力学。
b.在范围 \theta \in[-5.0,5.0] 内，哪个值 (\theta, \dot{\theta}) 代表不动点？哪些不动点是稳定的，哪些是不稳定的？
c.对于所有的 (\theta, \dot{\theta}) 对，我们的模型能够对动力学进行合理的估计吗？原因是什么？
