几种计算平面倒摆系统近似解的解析方法比较研究

Comparative Study on Several Analytical Methods for Approximate Solutions of Planar Inverted Pendulum System

摘要: 倒摆是常见的非线性动力系统，具有广泛的理论和实际研究价值，而采用解析方法计算倒摆的近似周期响应是常用研究手段。由于倒摆的运动过程中蕴含丰富的非线性动力学现象，如何选取合适的解析方法并进行特定周期行为的计算通常是能否进行有效讨论的关键。常见的定量解析分析方法包括参数扰动法、同伦分析法、能量法、平均法等。本文采用这四种方法分别对平面倒摆的控制方程进行求解，并采用Maple软件进行数值模拟对比，比较几种分析方法的优缺点及适用条件。本文对非线性倒摆系统进行的定量研究有助于理解倒摆的动力学行为，对类似的非线性动力系统的定量计算有较好的参考价值。

Abstract: The plane inverted pendulum system is a classical model with important theory value and vast application prospect. It is usual to research the periodic solutions of an inverted pendulum by using analytical methods. Since the inverted pendulum has strong nonlinearity and can exhibit very rich nonlinear phenomena, it is the key that how to choose an appropriate analytical method to calculate the periodic motions. There are four usual methods about calculating periodic solutions of a non-linear system: parameter perturbation method, homotopy analysis method, energy method and average method. In this paper, we will apply such four methods to analyze the periodic solutions of an inverted pendulum system, respectively. By comparing with numerical results obtained by using Maple software, we demonstrate the merits and demerits of the four methods as well as their application conditions. The research presented in this paper contributes to understanding nonlinear dynamics of an inverted pendulum and provides a good reference for calculating periodic solutions of other similar nonlinear systems.

1. 引言

倒摆是一个典型的形式简洁的非线性系统，其概念于20世纪60年代后期被正式提出，它在一定条件下可能会出现周期运动、旋转和混沌。倒摆具有巨大的研究价值，其运动的研究对于工程问题和科技发展有着很重大的意义，许多物理、工程、科研问题都可以简化为对平面倒摆模型的研究。如：人与结构的相互作用中如何保证人行走舒适性和稳定性的要求；双足机器人的直立行走过程中如何控制平衡；火箭发射或飞行过程中如何控制调整角度等。在研究倒摆的过程中，其近似周期响应一般采用解析方法计算。由于倒摆运动能反映出丰富的非线性动力学现象，如何选取合适的解析方法并进行特定周期行为的计算通常是能否进行有效讨论并解决问题的关键。

对倒摆系统的非线性振动问题的研究方法，主要分为实验法和分析法，其中分析方法又分为定性分析和定量分析。定性分析方法中常用的是1885年Poincaré H提出的相平面法，相平面法是一种可以用来研究非线性系统的解、平衡点及系统稳定性的图解方法。定量分析又可称作解析法，相较定性分析，其发展较为迅速，且种类也更加丰富。目前针对弱非线性问题的定量分析方法主要有：参数扰动法、多尺度法、坐标变换法等；针对强非线性问题的定量分析方法主要有：KBM法、能量法、平均法、伽辽金法、谐波平衡法、椭圆函数法、通论分析法、数值分析法等。2016年Ahsan等人 [1] 研究站立期间人体保持平衡稳定能力时建立了倒立单摆和双摆模型，考虑生理反映延迟的同时用伽勒金法分析模型，为校正器等辅助行走的器械的设计提供参考。同年徐等人 [2] 采用高斯白噪声扰动下的双连杆倒立摆模型研究人静止站立时如何维持平衡状态时，利用多尺度法分析其稳定性与大脑反应速度(时滞)的联系，指出反应迟钝的人更加容易失去静止站立的稳定性。2018年欧阳等人 [3] 选用扰动法求解基于平面倒摆模型的人行走动力学问题，阐述了行人保持稳定行走的内在机理。

本文首先建立平面倒摆物理模型，随后选用参数扰动法、同伦分析法、能量法、平均法这四种定量分析方法求解平面倒摆的运动方程，得到平面倒摆保持稳定运动的解，并对解析解和精确解的数值结果进行对比。最后比较这四种定量分析方法，阐述各自的优缺点。本文的分析有助于理解倒摆系统动力学行为，并为类似非线性动力系统的定量计算的方法选择提供参考。

2. 模型介绍

2.1. 物理模型

本文研究的平面倒立摆系统如图1所示，倒立摆的摆长、集中质量、摆与竖直方向的夹角和摆支点处附加的竖向周期作用分别为L、m、θ、。

Figure 1. The simple planer inverted pendulum model

图1. 简单平面倒立摆模型

2.2. 数学模型

如图1所示，该系统的动能和势能如下：

$T=\frac{1}{2}m\left[L\stackrel{\dot{}}{\theta }\cos {\left(\theta \right)}^2+{\left(\stackrel{\dot{}}{y}-L\stackrel{\dot{}}{\theta }\sin \left(\theta \right)\right)}^2\right],V=mg\left(y+L\cos \left(\theta \right)\right),$ (1)

由拉格朗日方程

(2)

式(2)中 $H=T-V$ ，将式(1)代入式(2)得到

$\ddot{\theta }+\left(\delta +\epsilon \cos \left(t\right)\right)\sin \left(\theta \right)=0,$ (3)

其中 $\epsilon =\frac{A}{L},\delta =-\frac{g}{L{\omega }^2}$ ，且 $0<\epsilon <1,-1<\delta <0$ 。本文考虑倒立摆摇摆角度大一些的情况，即 $\sin \left(\theta \right)\approx \theta -\frac{{\theta }^3}{6}$ ，在式(3)中用x代替θ，则有

$\ddot{x}+\left(\delta +\epsilon \cos \left(t\right)\right)\left(x-\frac{x^3}{6}\right)=0.$ (4)

3. 参数扰动法求解倒立摆系统

参数扰动法是一种常用的近似求解分析方法，它依赖于系统内参数或人工参数。参数扰动法的基本思路是，若扰动参数与系统的解无关，可以将解展开成该参数的幂级数的形式，从而将非线性问题转化为若干个线性子问题，然后逐次求解，每一步求解时得到的积分常数可由初始条件确定，最后取解的展开项的前几项来逼近系统的解析解。

利用参数扰动方法求解控制方程(4)，将δ展开成

$\delta ={\delta }_0+\epsilon {\delta }_1+{\epsilon }^2{\delta }_2+\cdots ,$ (5)

其中 ${\delta }_i$ 是待定常数，讨论解的稳定性和周期解需要讨论 ${\delta }_{\text{0}}=\text{0}$ 和 ${\delta }_0=\frac{1}{4}$ 这两种情况 [4] 。

当 ${\delta }_{\text{0}}=\text{0}$ 时，令解的形式为

$x=x_0+\epsilon x_1+{\epsilon }^2{x}_2+\cdots ,$ (6)

将式(5)，(6)代入控制方程式(4)，并使 ${\epsilon }^0$ ， ${\epsilon }^1$ ， ${\epsilon }^{\text{2}}$ 项的系数为0，则有

${\ddot{x}}_0+{\delta }_0\left(x_0-\frac{1}{6}{x}_0^3\right)=0,$ (7)

${\ddot{x}}_1+{\delta }_1\left(x_0-\frac{1}{6}{x}_0^3\right)+{\delta }_0{x}_1\left(1-\frac{1}{2}{x}_0^2\right)+\cos \left(t\right)\left(x_0-\frac{1}{6}{x}_0^3\right)=0,$ (8)

${\ddot{x}}_2-\frac{1}{2}{\delta }_1{x}_0^2{x}_1-\frac{1}{2}{\delta }_0{x}_0^2{x}_2-\frac{1}{2}{\delta }_0{x}_0{x}_1^2-\frac{1}{2}\cos \left(t\right)x_0^2{x}_1+\cos \left(t\right)x_1-\frac{1}{6}{\delta }_2{x}_0^3+{\delta }_2{x}_0+{\delta }_1{x}_1+{\delta }_0{x}_2=0.$ (9)

由于 ${\delta }_0=0$ ，(7)的解为 $x_0=a$ ，a是由初始条件来确定的待定常数。将 ${\delta }_0$ ， $x_0$ 的表达式代入式(8)得到

$x_1=-\frac{1}{6}\cos \left(t\right)a\left(a^2-6\right)+\frac{1}{2}\left(-a{\delta }_1+\frac{1}{6}{\delta }_1{a}^3\right)t^2+C_{1}t+C_2,$

为了消除久期项，令 $$ ，即 ${\delta }_1=0,C_1=C_2=0$ ，则

$x_1=-\frac{1}{6}\cos \left(t\right)a\left(a^2-6\right).$ (10)

将 $x_0,x_1,{\delta }_0,{\delta }_1$ 的表达式代入式(9)，并消除久期项，得到

$x_2=\frac{1}{96}a\left(12-8a^2+a^4\right)\left(\cos \left(2t\right)-1\right),$ (11)

${\delta }_2=\frac{1}{4}\left(a^2-2\right),$ (12)

因此，过原点的过渡曲线与其对应的周期解可以近似写为

$\delta =\frac{\text{1}}{\text{4}}\left(a^2-2\right){\epsilon }^{\text{2}},$ (13)

$x=a\left[1-\frac{1}{6}\left(a^2-6\right)\cos \left(t\right)\epsilon +\frac{1}{96}\left(12-8a^2+a^4\right)\left(\cos \left(2t\right)-1\right){\epsilon }^2\right].$ (14)

当 ${\delta }_0=\frac{1}{4}$ 时通过类似的计算过程可以近似得到延伸到 $\delta <0$ 的这一支分界线和其上对应的周期解

$\delta =\frac{1}{4}-\frac{1}{2}\epsilon +\left(\frac{1}{\text{32}}{a}^2-\frac{\text{1}}{\text{8}}\right){\epsilon }^2,$ (15)

$x=a\left\{\cos \left(\frac{t}{2}\right)\epsilon +\frac{1}{4}\cos \left(\frac{3t}{2}\right){\epsilon }^2+\left(\frac{1}{48}\cos \left(\frac{5t}{2}\right)-\left(\frac{1}{192}{a}^2+\frac{1}{16}\right)\cos \left(\frac{3t}{2}\right){\epsilon }^3\right)\right\}.$ (16)

分别取ε = 0.1、0.2、0.5对系统(4)进行数值模拟。考虑到初值对(4)的结果有影响，在模拟计算时对每个ε分别取阿a = 0.05、0.5、1。数值结果和解析结果(16)的比较如图2所示。

$\epsilon =0.1a=0.05$ $\epsilon =0.1a=0.5$ $\epsilon =0.1a=1$ $\epsilon =0.2,a=0.05$ $\epsilon =0.2,a=0.5$ $\epsilon =0.2,a=1$ $\epsilon =0.5,a=0.05$ $\epsilon =0.5,a=0.5$ $\epsilon =0.5,a=1$

Figure 2. The comparison between the numerical results of Equation (4) and the analytical results of Equation (24), solid line is the numerical results, dot line is the analytical results

图2. 方程(4)的数值解和解析解(24)结果比较，实线为数值结果，点代表解析解

显然，a和ε较大时都会降低近似解析解的精度。如果要提高精度需要计算ε的更高阶项。图3中给出了 $\delta <0$ 时初始条件a对分界线(13)和(15)的影响。随着初始摆角的增加，倒摆的稳定性区域不断减少。从参数平面上δ-轴上 $\delta =0$ 点出发的分界线更易受到初始条件的影响。

Figure 3. The influence of the initial condition a on the transition curves (13) and (15)

图3. 初始条件a对参数平面上分界线(13)和(15)的影响

4. 同伦分析法求解倒立摆系统

“同伦”作为拓扑学中的一个重要的概念，是于1911年由brouwer首次提出。1992年，廖世俊开创性地利用同伦分析法来处理非线性问题，并创立了“同伦分析方法”，随后几十年中廖继续完善同伦分析法，改进后的同伦分析法也已被验证满足近似方法的评判标准，目前存在的一个公认的评价解析解近似方法优劣的标准如下：

1) 是否能将一个非线性问题转化为诸多线性子问题；

2) 是否能使线性子问题快速地被求解；

3) 是否能让线性子问题的解组成的近似解是逼近原始问题的，即所得解析结果是否足够精确。

在同伦分析法中，我们可以这样表述，设 $f,g:X\to Y$ 都是连续映射， $I=\left[0,1\right]$ ，若存在连续映射 $H:X\times I\to Y$ ，使得所有的 $x\in X$ 时有， $H\left(x,0\right)=f\left(x\right)$ ，，则称H为从f到g的一个同伦，此时 $f_0=f$ ， $f_{\text{1}}=g$ 。

基于同伦的基本思想，即数学上的连续变形，我们可以试采用同伦分析法 [5] 解决强非线性问题。在控制方程(4)中，若用Ω表示最后求解的周期解的频率，令 $\tau =\Omega t$ ，则有

$\gamma \ddot{x}+\left(\delta +\epsilon \cos \left(\tau \right)\right)\left(x-\frac{1}{6}{x}^3\right)=0,$ (17)

其中“·”表示对 $\tau$ 求导， $\gamma ={\Omega }^2$ 。

由同伦分析方法，构造下述方程组，

$\left(1-q\right)L\left[x\left(\tau ;q\right)-x_0\left(\tau \right)\right]=hqN\left[x\left(\tau ;q\right)\right],$ (18)

其中嵌入变量 $q\in \left[0,1\right]$ ，h为非零辅助参数， $x_0\left(\tau \right)=x\left(\tau ;0\right)$ 是解的初始猜测解，L是线形辅助算子，N是非线性算子。基于(17)，非线性算子 $$ 可以写作

$N\left[x\left(\tau ;q\right)\right]=\gamma \left(q\right)\frac{{\partial }^2}{\partial {\tau }^2}x\left(\tau ;q\right)+\left(\delta +\epsilon \cos \left(\tau \right)\right)\left(x\left(\tau ;q\right)-\frac{1}{6}x{\left(\tau ;q\right)}^3\right),$ (19)

当q在 $\left[0,1\right]$ 间变化时， $\gamma \left(q\right)$ 也被定义为一个同伦。将(19)代入(18)得到

$\begin{array}{l}\left(1-q\right)L\left[x\left(\tau ;q\right)-x_0\left(\tau \right)\right]\\ =hq\left[\gamma \left(q\right)\frac{{\partial }^2}{\partial {\tau }^2}x\left(\tau ;q\right)+\left(\delta +\epsilon \cos \left(\tau \right)\right)\left(x\left(\tau ;q\right)-\frac{1}{6}x{\left(\tau ;q\right)}^3\right)\right],\end{array}$ (20)

这就是零阶形变方程，其对应的初始条件为

$x\left(0;q\right)=x_0\left(0\right)=c,\stackrel{\dot{}}{x}\left(0;q\right)=0.$ (21)

将 $x\left(\tau ;q\right)$ 和 $\gamma \left(q\right)$ 分别展开成

$x\left(\tau ;q\right)=x_0\left(\tau \right)+{\displaystyle \underset{n=1}{\overset{+\infty }{\sum }}{x}_n}\left(\tau \right)q^n,$ (22)

$x_n\left(\tau \right)={\frac{1}{n!}\frac{{\partial }^n}{\partial q^n}x\left(\tau ;q\right)|}_{q=0},$

$\gamma \left(q\right)={\gamma }_0+{\displaystyle \underset{n=1}{\overset{+\infty }{\sum }}{\gamma }_n}\left(\tau \right)q^n,$ (23)

${\gamma }_n={\frac{1}{n!}\frac{{\text{d}}^n}{\text{d}{q}^n}\gamma \left(q\right)|}_{q=0}.$

先将零阶形变方程式(20)的左右两侧同时对q求m阶导数，然后等式两侧同时除以 $m!$ ，令 $$ ，则可得到m阶形变方程

$\begin{array}{l}L\left[\frac{1}{m!}\frac{{\partial }^m}{\partial q^m}x\left(\tau ;q\right)\right]-L\left[\frac{1}{\left(m-1\right)!}\frac{{\partial }^{m-1}}{\partial q^{m-1}}x\left(\tau ;q\right)\right]\\ =\frac{h}{\left(m-1\right)!}{D}^{m-1}N\left[x\left(\tau ;q\right)\right],\end{array}$ (24)

式(24)又可表述成

$L\left[x_m\left(\tau \right)-{\chi }_m{x}_{m-1}\left(\tau \right)\right]=h{R}_{m}N\left[x_{m-1},{\gamma }_{m-1}\right],$ (25)

其中

采用同伦分析方法求解问题时，有很大的自由去选择基函数和线形辅助算子，令基函数为如下形式

(26)

按照基函数的形式，初始猜测解可以写作 $x_0\left(\tau \right)=c\cos \left(\tau \right)$ ，考虑到方程(17)的二阶性质，线性辅助算子 $L\left[x\left(\tau ;q\right)\right]$ 可以表述为

$L\left[x\left(\tau ;q\right)\right]=\frac{{\partial }^2}{\partial {\tau }^2}x\left(\tau ;q\right)+x\left(\tau ;q\right),$ (27)

具有性质

$L\left[C_1\sin \left(\tau \right)+C_2\cos \left(\tau \right)\right]=0.$

结合式(25)、(27)，得到

$$ (28)

当 $m=1,2,\cdots$ 时，由(28)计算得到

$\begin{array}{l}{x}_1\left(\tau \right)=\frac{1}{6}hc(c^2\epsilon \cos {\left(\tau \right)}^5+c^2\left(\delta -\epsilon \right)\cos {\left(\tau \right)}^4-c^2\delta \cos {\left(\tau \right)}^3+\left(6-\frac{3}{4}{c}^2\right)\delta \cos {\left(\tau \right)}^2\\ +\left(\frac{3}{4}{c}^2\delta +3\tau \epsilon \sin \left(\tau \right)-6\delta \right)\cos \left(\tau \right)+3\delta \tau \sin \left(\tau \right)),\\ x_{\text{2}}\left(\tau \right)=f_2\left(\tau \right),\\ \cdots \end{array}$ (29)

$\begin{array}{l}{\gamma }_0=-\frac{1}{8}\delta \left(c^2-8\right),\\ {\gamma }_1=\frac{1}{64}{\delta }^2\left(c^2-8\right)\left(9c^2-8\right)h,\\ {\gamma }_2=\cdots ,\end{array}$ (30)

其中， $x_2\left(\tau \right)$ 的具体表达式过于冗长，这里直接用表示。根据式(22)、(23)，令 $$ ，可以得到解析解的形式为

$\begin{array}{l}x\left(\tau \right)=x_0\left(\tau \right)+x_1\left(\tau \right)+x_2\left(\tau \right)+\cdots ,\\ \gamma ={\gamma }_0+{\gamma }_1+\cdots .\end{array}$ (31)

将上述得到的、 $x_1\left(\tau \right)$ 、 $x_{\text{2}}\left(\tau \right)$ 的表达式和式(30)代入式(31)中，可以得到原始控制方程的解：

$\begin{array}{l}x\left(\tau \right)=x_0\left(\tau \right)+x_1\left(\tau \right)+x_{\text{2}}\left(\tau \right)+\cdots \\ =c\cos \left(\tau \right)+[\frac{1}{6}hc(c^2\epsilon \cos {\left(\tau \right)}^5+c^2\left(\delta -\epsilon \right)\cos {\left(\tau \right)}^4-c^2\delta \cos {\left(\tau \right)}^3+\left(6-\frac{3}{4}{c}^2\right)\delta \cos {\left(\tau \right)}^2\\ +\left(\frac{3}{4}{c}^2\delta +3\tau \epsilon \sin \left(\tau \right)-6\delta \right)\cos \left(\tau \right)+3\delta \tau \sin \left(\tau \right))]+f_{\text{2}}\left(\tau \right)+\cdots ,\\ \gamma ={\gamma }_0+{\gamma }_1+\cdots =-\frac{1}{8}\delta \left(c^2-8\right)+\frac{1}{64}{\delta }^2\left(c^2-8\right)\left(9c^2-8\right)h+\cdots .\end{array}$ (32)

由式(32)可知，解的表达式中包含多个参数，其中h是非零的收敛控制参数，用来控制结果的收敛性，c是初始位移，ε和δ之间没有很明确的联系，下面将通过数值模拟给出参数c和ε对于解的影响。当 $\delta =-0.2,h=-1/1000,c=0.005$ 时分别选取 $\epsilon =\text{0}\text{.01},0.05,0.08$ 对(32)中第一式的解析解和对应于(17)的精确解进行数值比较，其结果如图4所示。随着ε的增大，同伦法的精度迅速下降。

当 $\epsilon =\text{0}\text{.01},\delta =-0.2,h=-1/1000$ 时，选取 $$ 进行比较，如图5所示。随着c的增大，拟合效果并没有很明显的变化。

(1) $\epsilon =0.01$ (2) $\epsilon =0.0\text{5}$ (3) $\epsilon =0.0\text{8}$

Figure 4. The comparison between analytical and numerical results of Equation (32) when $\delta =-0.2$ , $h=-1/1000$ , $c=0.005$ , solid line is the numerical results, dot line is the analytical results

图4. 解析解(32)和数值解的比较，其中 $\delta =-0.2$ ， $h=-1/1000$ ， $c=0.005$ 实线为数值解，点代表解析解

(1) $c=0.005$ (2) $c=0.05$ (3) $c=0.5$

Figure 5. The influence of c on the analytical and numerical results of Equation (32) when $\epsilon =\text{0}\text{.01,}\delta =-0.2,h=-1/1000$ , solid line is the numerical results, dot line is the analytical results

图5. c改变对解析解(32)的影响，其中 $\epsilon =\text{0}\text{.01,}\delta =-0.2,h=-1/1000$ ，实线为数值解，点代表解析解

5. 能量法求解倒立摆系统

能量法能很好地处理强非线性问题从而得到周期解 [6] 。基于能量法的基本思想，可以知道如果物体做周期运动，在一个周期的时间里对能量进行平均，得到的值应该是一个定常数。如果物体的周期运动是渐近稳定的，则在领域中，即与周期相同的时间内，对运动体的能量进行平均得到的数值趋近于周期运动的平均能量。

考虑到控制方程(4)是强非线性非自治系统，将其中仅含x的所有项记作 $g\left(x\right)$ ，其它各项记作 $f\left(x,t\right)$ ，则有关于方程(4)的正则形式为

$\ddot{x}+g\left(x\right)+f\left(x,t\right)=0,$ (33)

其中

$g\left(x\right)=\delta \left(x-\frac{x^3}{6}\right),$ (34)

$f\left(x,t\right)=\epsilon \cos \left(t\right)\left(x-\frac{x^3}{6}\right).$ (35)

注意到连续， $g\left(x\right)$ 连续，且满足条件 $xg{\left(x\right)}_{x\ne 0}>0$ 。 $g\left(x\right)$ 对应的势能及系统(33)对应的总能量分别为

$V\left(x\right)={\displaystyle {\int }_0^{x}g\left(\xi \right)\text{d}\xi }=\frac{\delta }{2}{x}^2-\frac{\delta }{24}{x}^4,$ (36)

$E\left(x,\stackrel{\dot{}}{x}\right)=\frac{1}{2}{\stackrel{\dot{}}{x}}^2+V\left(x\right)=\frac{1}{2}{\stackrel{\dot{}}{x}}^2+\frac{\delta }{2}{x}^2-\frac{\delta }{24}{x}^4.$ (37)

根据能量法，系统(43)对应的能量坐标变换公式如下：

$x=a\cos \left(\theta \right),$ (38)

$\begin{array}{c}\stackrel{\dot{}}{x}=-\sqrt{2\left(V\left(a\right)-V\left(a\cos \left(\theta \right)\right)\right)}\\ =-a\sin \left(\theta \right)\sqrt{\delta }\left(1-\frac{a^2}{16}-\frac{a^2}{48}\cos \left(2\theta \right)\right).\end{array}$ (39)

将方程(33)的各项乘以 $\stackrel{\dot{}}{x}$ ，再结合(37)可以得到

$\begin{array}{l}\stackrel{\dot{}}{x}\ddot{x}+\stackrel{\dot{}}{x}g\left(x\right)=-\stackrel{\dot{}}{x}f\left(x,t\right),\\ \stackrel{\dot{}}{x}\ddot{x}+\stackrel{\dot{}}{x}g\left(x\right)=\frac{\text{d}}{\text{d}t}\left(\frac{1}{2}{\stackrel{\dot{}}{x}}^2+V\left(x\right)\right)=\frac{\text{d}E}{\text{d}t},\end{array}$ (40)

则

$\frac{\text{d}E}{\text{d}t}=-\stackrel{\dot{}}{x}f\left(x,t\right),$ (41)

将式(35)，(38)，(39)代入式(41)，整理可以得到

$\frac{\text{d}E}{\text{d}t}=\frac{a^2}{4}\sqrt{\delta }\epsilon \left(1-\frac{a^2\cos {\left(\theta \right)}^2}{6}\right)\left(\sin \left(2\theta +t\right)+\sin \left(2\theta -t\right)\right)\left(1-\frac{a^2}{16}-\frac{a^2}{48}\cos \left(2\theta \right)\right)=F\left(a,\theta ,t\right),$ (42)

这就是系统(33)对应的能量坐标系中的第一个方程。

将式(38)对t求导，结合式(39)，加以考虑a是E的函数，得到

$\frac{\text{d}\theta }{\text{d}t}=\frac{1}{a\sin \left(\theta \right)}\left[\frac{\text{d}a}{\text{d}E}\cos \left(\theta \right)\frac{\text{d}E}{\text{d}t}+\sqrt{2\left(V\left(a\right)-V\left(a\cos \left(\theta \right)\right)\right)}\right].$ (43)

由式(36)可知 $V\left(x\right)$ 为偶函数，则有

$\begin{array}{l}V\left(-a+b\right)=E,\\ V\left(a+b\right)=E,\end{array}$

将这两个式子对E求导，则可以求出 $$ 和 $\frac{\text{d}b}{\text{d}E}$ ，这里令 $b=0$ ，则

$\frac{\text{d}a}{\text{d}E}=\frac{g\left(-a\right)-g\left(a\right)}{2g\left(-a\right)g\left(a\right)},$ (44)

将式(42)、(44)代入式(43)，经计算得到

$\frac{\text{d}\theta }{\text{d}t}=\sqrt{\delta }\left(1-\frac{a^2}{16}-\frac{a^2}{48}\cos \left(2\theta \right)\right)\left(1+\frac{\epsilon \cos {\left(\theta \right)}^2\cos \left(t\right)\left(1-\frac{a^2\cos {\left(\theta \right)}^2}{6}\right)}{\delta \left(1-\frac{a^2}{6}\right)}\right)\ne \text{0},$ (45)

这就是系统(33)对应的能量坐标系中的第二个方程。

直接考虑，则 $t=2\left(\theta -\varphi \right)$ ，将其代入式(42)和(45)，得到

$\{\begin{array}{l}\frac{\text{d}E}{\text{d}t}=\frac{a^2}{4}\sqrt{\delta }\epsilon \left(1-\frac{a^2\cos {\left(\theta \right)}^2}{6}\right)\left(\sin \left(\text{4}\theta -2\varphi \right)+\sin \left(2\varphi \right)\right)\left(1-\frac{a^2}{16}-\frac{a^2}{48}\cos \left(2\theta \right)\right),\\ \frac{\text{d}\theta }{\text{d}t}=\sqrt{\delta }\left(1-\frac{a^2}{16}-\frac{a^2}{48}\cos \left(2\theta \right)\right)\left(1+\frac{\left(\cos \left(2\varphi \right)+2\cos \left(\text{2}\theta -2\varphi \right)+\cos \left(\text{4}\theta -2\varphi \right)\right)\frac{\epsilon }{4}\left(1-\frac{a^2\left(\cos \left(2\theta \right)+1\right)}{12}\right)}{\delta \left(1-\frac{a^2}{6}\right)}\right)\end{array},$ (46)

结合上式，要使周期解存在，则有

$\{\begin{array}{l}\frac{a^2}{4}\sqrt{\delta }\epsilon \left(1-\frac{a^2}{\text{12}}\right)\sin \left(2\varphi \right)\left(1-\frac{a^2}{16}\right)=\text{0},\\ \sqrt{\delta }\left(1-\frac{a^2}{16}\right)+\left(1-\frac{a^2}{16}\right)\frac{\epsilon \cos \left(2\varphi \right)\left(1-\frac{a^2}{12}\right)}{4\sqrt{\delta }\left(1-\frac{a^2}{6}\right)}=\frac{1}{2},\end{array}$ (47)

求解(47)第一式，得到 $\varphi =0$ ， $\varphi =\frac{\text{π}}{\text{2}}$ ，将其代入(47)第二式，则有

$\sqrt{\delta }\left(1-\frac{a^2}{16}\right)\left(1+\frac{\epsilon \left(1-\frac{a^2}{12}\right)}{4\delta \left(1-\frac{a^2}{6}\right)}\right)=\frac{1}{2},\varphi =\text{0}$ (48)

$\sqrt{\delta }\left(1-\frac{a^2}{16}\right)\left(1-\frac{\epsilon \left(1-\frac{a^2}{12}\right)}{4\delta \left(1-\frac{a^2}{6}\right)}\right)=\frac{1}{2},\varphi =\frac{\text{π}}{\text{2}}$ (49)

将 $\varphi =0$ 和 $\varphi =\frac{\text{π}}{\text{2}}$ 分别代入式(46)得到

$\{\begin{array}{l}\frac{\text{d}E}{\text{d}t}=\frac{a^2}{4}\sqrt{\delta }\epsilon \left(1-\frac{a^2\cos {\left(\theta \right)}^2}{6}\right)\sin \left(\text{4}\theta \right)\left(1-\frac{a^2}{16}-\frac{a^2}{48}\cos \left(2\theta \right)\right),\\ \frac{\text{d}\theta }{\text{d}t}=\sqrt{\delta }\left(1-\frac{a^2}{16}-\frac{a^2}{48}\cos \left(2\theta \right)\right)\left(1+\frac{\left(\text{1}+2\cos \left(\text{2}\theta \right)+\cos \left(\text{4}\theta \right)\right)\frac{\epsilon }{4}\left(1-\frac{a^2\left(\cos \left(2\theta \right)+1\right)}{12}\right)}{\delta \left(1-\frac{a^2}{6}\right)}\right),\end{array}$ (50)

和

$\{\begin{array}{l}\frac{\text{d}E}{\text{d}t}=\frac{a^2}{4}\sqrt{\delta }\epsilon \left(1-\frac{a^2\cos {\left(\theta \right)}^2}{6}\right)\left(-\sin \left(\text{4}\theta \right)\right)\left(1-\frac{a^2}{16}-\frac{a^2}{48}\cos \left(2\theta \right)\right),\\ \frac{\text{d}\theta }{\text{d}t}=\sqrt{\delta }\left(1-\frac{a^2}{16}-\frac{a^2}{48}\cos \left(2\theta \right)\right)\left(1+\frac{\left(-1-2\cos \left(\text{2}\theta \right)-\cos \left(\text{4}\theta \right)\right)\frac{\epsilon }{4}\left(1-\frac{a^2\left(\cos \left(2\theta \right)+1\right)}{12}\right)}{\delta \left(1-\frac{a^2}{6}\right)}\right).\end{array}$ (51)

式(38)两边求导可以写成

$\stackrel{\dot{}}{x}=\frac{\text{d}a}{\text{d}E}\frac{\text{d}E}{\text{d}t}\cos \left(\theta \right)-a\sin \left(\theta \right)\frac{\text{d}\theta }{\text{d}t},$ (52)

将式(44)，式(50)或(51)代入式(52)之后，进行积分运算即可得到式(33)在 $\varphi =0$ 或 $\varphi =\frac{\text{π}}{\text{2}}$ 时的解。

分别在 $\varphi =0$ 或 $\varphi =\frac{\text{π}}{\text{2}}$ 的情况下，取 $\epsilon =0.00\text{1,0}\text{.01,0}\text{.1}$ ，a = 0.5(a可视为振幅)对(4)进行数值计算并和解析解比较，结果如图6所示。

(1) $\varphi =0,\epsilon =0.001$ (2) $\varphi =0,\epsilon =0.01$ (3) $\varphi =0,\epsilon =0.1$

Figure 6. The comparison between analytical and numerical results with energy method when $a=\text{0}.5$ , solid line is the numerical results, dot line is the analytical results

图6.时能量法解析解和数值解比较，实线为数值解，点代表解析解

(1) $\varphi =0,a=0.01$ (2) $\varphi =0,a=0.1$ (3) $\varphi =0,a=0.4$ (4) $\varphi =\frac{\text{π}}{2},a=0.01$ (5) $\varphi =\frac{\text{π}}{2},a=0.1$ (6) $\varphi =\frac{\text{π}}{2},a=0.4$

Figure 7. The comparison between analytical and numerical results with energy method when $\epsilon =\text{0}\text{.1}$ , solid line is the numerical results, dot line is the analytical results

图7. $\epsilon =\text{0}\text{.1}$ 时能量法解析解和数值解比较，实线为数值解，点代表解析解

另外，分别在 $\varphi =0$ 或 $\varphi =\frac{\text{π}}{\text{2}}$ 的情况下，取 $a=\text{0}\text{.01,0}\text{.1,0}\text{.4}$ ， $\epsilon =0.1$ 进行数值计算并与解析解比较，结果如图7所示。随着a增大，解析解与数值解的差异逐渐增大，这和扰动法的变化趋势是一样的。

6. 平均法求解倒立摆系统

平均法是一种简便的定量分析方法，在平均法中，我们首先考虑常数变易法，将系统解的振幅和相位均表示成时间 $t$ 的函数，再加以利用慢变参数法 [7] ，即系统解的振幅和相位均随着时间缓慢变化，则可以把问题转变为含有振幅和相位的关系式。分析关于振幅和相位的关系式，就可以得到系统周期解时对应的参数关系曲线。

采用平均法 [8] 求解控制方程(4)时，其解可以写成如下形式

$x\left(t\right)=B\left(t\right)\cos \left(\theta \left(t\right)\right)=B\left(t\right)\cos \left(\frac{1}{2}-\varphi \left(t\right)\right),$ (53)

令

$\stackrel{\dot{}}{B}\cos \left(\theta \right)+B\stackrel{\dot{}}{\varphi }\sin \left(\theta \right)=0,$ (54)

将式(53)、(54)代入方程(4)得到

$-\frac{1}{2}B\cos \left(\theta \right)\stackrel{\dot{}}{\theta }-\frac{1}{2}\stackrel{\dot{}}{B}\sin \left(\theta \right)+\delta B\cos \left(\theta \right)+f\left(B\cos \left(\theta \right),t\right)=0,$ (55)

其中

$f\left(B\cos \left(\theta \right),t\right)=\epsilon B\cos \left(t\right)\cos \left(\theta \right)-\frac{1}{6}{\left(B\cos \left(\theta \right)\right)}^3\left(\delta +\epsilon \cos \left(t\right)\right).$

将式(54)乘以 $\sin \left(\theta \right)$ ，式(55)乘以 $\cos \left(\theta \right)$ ，并把二者相加，可以得到

$\left(\delta -\frac{1}{4}\right)B\cos {\left(\theta \right)}^2+\frac{1}{2}B\stackrel{\dot{}}{\varphi }+\cos \left(\theta \right)f\left(B\cos \left(\theta \right),t\right)=0,$ (56)

将式(55)乘以 $\sin \left(\theta \right)$ ，式(54)乘以 $\cos \left(\theta \right)$ ，用前者减去后者，可以得到

$\left(\delta -\frac{1}{4}\right)B\cos \left(\theta \right)\sin \left(\theta \right)-\frac{1}{2}\stackrel{\dot{}}{B}+\sin \left(\theta \right)f\left(B\cos \left(\theta \right),t\right)=0.$ (57)

若B、 $\varphi$ 均随着时间缓慢变化，那么可以利用慢变参数法 [8] ，即在一个周期里B、 $\varphi$ 可以当作定常数。将式(56)、(57)分别对θ在0到 $\text{2π}$ 上进行积分，则有

$B\left(\stackrel{\dot{}}{\varphi }+\delta -\frac{1}{4}-\frac{1}{8}\delta B^2+\frac{1}{2}\epsilon \cos \left(2\varphi \right)-\frac{1}{12}\epsilon B^2\cos \left(2\varphi \right)\right)=0,$ (58)

(59)

令 $\stackrel{\dot{}}{B}=0$ 、 $\stackrel{\dot{}}{\varphi }=0$ ，使得求解出来的运动形式为稳定的振动，考虑式(58)和(59)，得到

1) $B=0$ ；

2) $B\ne 0$ ， $2\varphi =0$ ， $\delta -\frac{1}{4}-\frac{1}{8}\delta B^2+\frac{1}{2}\epsilon -\frac{1}{12}\epsilon B^2=0$ ；

3) $B\ne 0$ ， $2\varphi =\text{π}$ ， $\delta -\frac{1}{4}-\frac{1}{8}\delta B^2-\frac{1}{2}\epsilon +\frac{1}{12}\epsilon B^2=0$ ；

Figure 8. The influence of B on the parameter curves when $\varphi =\text{0}$

图8. $\varphi =\text{0}$ 时B对参数曲线的影响

这里考虑B为正数的情况。因为 $\delta <0$ ，所以分析上述第二种情况，即 $2\varphi =0$ 。如图8所示，当 $2\varphi =0$ 时，分别取、0.5、0.8、1，随着B的增大，参数曲线上移。在图8中的4条参数曲线上标记16个点， $P_{ai}$ 、 $P_{bi}$ 、 $$ 、 $P_{di}$ ( $i=1,2,3,4$ )，从中选取9个点进行数值模拟，结果如图9所示。B对解析解精度的影响很小，ε的大小对解析解的精度起到关键性的作用，ε越小解析解和数值结果的吻合程度越好。

$\epsilon =0.1,P_{b1},B=0.05$ $\epsilon =0.1,P_{b1},B=0.5$ $\epsilon =\text{0}\text{.1},P_{b4},B=1$ $\epsilon =\text{0}\text{.2},P_{c1},B=0.05$ $\epsilon =\text{0}\text{.2},P_{c2},B=0.5$ $\epsilon =\text{0}\text{.2},P_{c4},B=1$ $\epsilon =\text{0}\text{.3},P_{d1},B=0.05$ $\epsilon =\text{0}\text{.3},P_{d2},B=0.5$ $\epsilon =\text{0}\text{.3},P_{d2},B=1$

Figure 9. The comparison between analytical and numerical results with averaging method when ε = 0.1,0.2,0.3 , solid line is the numerical results, dot line is the analytical results

图9. 当 $\epsilon =\text{0}\text{.1},\text{0}\text{.2},\text{0}\text{.3}$ 平均法解析解和数值解比较，实线为数值解，点代表解析解

为了方便比较四种解析方法求解倒立摆稳定周期解的效果，取初始位移 $x\left(0\right)=0.5$ ， $\epsilon =0.1,\text{0}\text{.3}$ 分别绘制四种解析方法的结果曲线图，具体结果见图10。可以看出在ε较小时四种方法得到的解析解表达式精度相差不大。当 $\epsilon >0.1$ 时参数扰动法和平均法比同伦法和能量法的精度高。

参数扰动法 能量法 平均法 同伦分析法 (1) $\epsilon =0.1,x\left(0\right)=0.5$ 参数扰动法 能量法 平均法 同伦分析法(2) $\epsilon =0.3,x\left(0\right)=0.5$

Figure 10. The comparison between analytical and numerical results with the four analysis methods: parameter perturbation method, energy method, averaging method and homotopy analysis method, when $\epsilon =\text{0}\text{.1},\text{0}\text{.2},x\left(0\right)=0.5$ , solid line is the numerical results, dot line is the analytical results

图10. 四种分析方法的解析解和数值解比较结果， $\epsilon =\text{0}\text{.1},\text{0}\text{.2},x\left(0\right)=0.5$ ，实线为数值解，点代表解析解

7. 结论

本文采用四种解析方法定量计算了平面倒立摆系统的稳定周期解，并将解析解和数值解比较展示了解析解的精度及影响精度的因素。对比四种方法有如下结论：

1) 参数扰动法、同伦分析法、能量法和平均法都能近似得到平面倒摆稳定周期解的解析表达式。表达式的精度受初始条件和系统控制方程中非线性项系数大小的影响。当摆的初始摆角增加或者非线性项系数增加时解析解的精度下降较快。

2) 采用参数扰动法、能量法、平均法研究倒摆能直接得到周期解出现的参数条件，而同伦法不能。采用同伦分析方法时，基函数、初始猜测解和线性辅助算子的选择主观性较强，具有自由灵活的特征，但这也为求解周期解带来确定参数的困难。

3) 在 $\text{0}<\epsilon <0.1$ 范围内四种方法得到的解析解表达式精度相差不大，其中平均法求解过程最为简洁；当 $\epsilon >0.1$ 时参数扰动法和平均法比同伦法和能量法的精度高。如果要提高精度需要计算关于ε的高次项。四种方法中采用参数扰动法计算ε的更高次项相对而言计算量最小，能量法、平均法计算较为复杂。同伦法计算过程复杂的同时还存在参数难以确定的困难。

基金项目

本文感谢国家自然科学基金(11472160)的支持。