直线二阶倒立摆之数学建模

前言

本文是关于二阶倒立摆模型如何在平衡点处得到其的线性化模型的说明，介绍的比较详细，但是大部分内容参考于前人所写的论文，这些参考文献，我将列在本文的最后，但是我所建立的模型不同于论文中的模型，基本上计算过程我也是手算了一遍，在这里写出来，也是为了给想要计算一遍的人一个参考。如果大家迫切需要得到结论，大家可以直接照抄论文内的模型，论文内的模型更具有普适性。我这个模型是个特例。

仿真物理模型的建立

为了达到更加逼真以及生动的演示效果，我将采用simulink中的simscape模块对二阶直线倒立摆进行建模。关于如何使用simscape对二阶倒立摆建模我将在另外一篇文章中进行讲解。直接给出建立的模型。如下

物理意义描述为

物理含义	符号
滑轨质量	M
滑轨位置	x
摆杆1质量	m1
摆杆1角度	θ1
摆杆1半杆长	l1
摆杆2质量	m2
摆杆2角度	θ2
摆杆2半杆长	l2

真正完整的真实的倒立摆模型在两个摆杆之间的应当还会有一个转轴，但是我这里是使用simulink仿真，并不是为了现实的控制，所以直接忽略了，但是在后续会给出完整模型的参数

二阶倒立摆的数学模型

使用拉格朗日方程建模

使用拉格朗日方程对二阶倒立摆进行数学模型的建模。拉格朗日方程为
$$\frac{d}{dt}\frac{\partial L}{\partial {\dot{q}}_i}-\frac{\partial L}{\partial q_i}=f_i(i=1,2,3...)\text{\hspace{1em}(1)}$$

此处关于拉格朗日方程的详细说明大家可以参照小时物理百科——一个页面非常精美的物理概念的集合网站，拉格朗日概念传送门

其中q_i为系统的广义坐标，f_i为系统第i个广义坐标上的系统外力，L为系统的拉格朗日量，在二阶倒立摆中，描述系统的广义坐标即为(x,θ₁,θ₂)；而显然在θ₁和θ₂上系统的外力为零；而拉格朗日量表述为：
$$L(q_1,q_2,...,{\dot{q}}_1,{\dot{q}}_2,...,t)=T(q_1,q_2,...,{\dot{q}}_1,{\dot{q}}_2,...,t)-V(q_1,q_2,...,t)\text{\hspace{1em}(2)}$$
其中T与V分别为系统的动能与势能

数学模型的建立

在二阶倒立摆中，系统动能表述为：
$$T=T_0+T_1+T_2\text{\hspace{1em}(3)}$$
其中T₀、T₁、T₂分别为小车、摆杆1、摆杆2的动能，下面我们再计算T₀、T₁、T₂
$$T_0=\frac{1}{2}M{\dot{x}}^2\text{\hspace{1em}(4)}$$
$$T_1=\frac{1}{6}{m}_1{l}_1^2{\dot{\theta }}_1^2+\frac{1}{2}{m}_1\left[{(\frac{d}{dt}\left(x-l_{1}sin{\theta }_1\right))}^2+{\left(\frac{d}{dt}{l}_{1}cos\dot{{\theta }_1}\right)}^2\right]\text{\hspace{1em}(5)}$$
$$T_2=\underset{转动动能}{\underbrace{\frac{1}{6}{m}_2{l}_2^2{\dot{\theta }}_2^2}}+\frac{1}{2}{m}_2\left[\underset{x方向位置的导数，即x方向的速度}{\underbrace{{\left(\frac{d}{dt}\left(x-2l_{1}sin{\theta }_1-l_{2}sin{\theta }_2\right)\right)}^2}}+\underset{y方向位置的导数，即y方向的速度}{\underbrace{{\left(\frac{d}{dt}\left(2l_{1}cos{\theta }_1+l_{2}cos{\theta }_2\right)\right)}^2}}\right]\text{\hspace{1em}(6)}$$

这里简单的解释以下，T₁、T₂如何求的。即为转动动能和平动动能，平动动能里再细分为x和y方向的动能，上面这个公式渲染的不是很好，看来csdn对latex的支持也不太好

而系统势能描述为
$$V=V_0+V_1+V_2\text{\hspace{1em}(7)}$$
同样为小车、摆杆1、摆杆2的势能之和
$$V_0=0\text{\hspace{1em}(8)}$$
$$V_1=m_{1}g{l}_{1}cos{\theta }_1\text{\hspace{1em}(9)}$$
$$V_2=m_{2}g\left(2l_{1}cos{\theta }_1+l_{2}cos{\theta }_2\right)\text{\hspace{1em}(10)}$$
简单整合一下
$$\begin{array}{rl}T=& \frac{1}{2}\left(M+m_1+m_2\right){\dot{x}}^2+\left(\frac{2}{3}{m}_1+2m_2\right)l_1^2{\dot{\theta }}_1^2\\ & +\frac{2}{3}{m}_2{l}_2^2{\dot{\theta }}_2^2-\left(m_1+2m_2\right)l_1\dot{x}{\dot{\theta }}_{1}cos{\theta }_1\\ & -m_2{l}_2\dot{x}{\dot{\theta }}_{2}cos{\theta }_2+2m_2{l}_1{l}_2{\dot{\theta }}_1{\dot{\theta }}_{2}cos({\theta }_1-{\theta }_2)\end{array}\text{\hspace{1em}(11)}$$
$$V=m_{1}g{l}_{1}cos{\theta }_1+m_{2}g(2l_{1}cos{\theta }_1+l_{2}cos{\theta }_2)\text{\hspace{1em}(12)}$$

这里简单的说明一下，大家如果在建模的时候所选取的角度的正负不一样的话，这里有一些符号会又变化，本公式是以摆杆在铅垂线的左侧时的角度为正的。

前面已经提到，拉格朗日算子为
$$L=T-V\text{\hspace{1em}(13)}$$
而在该二阶摆系统中，在广义坐标θ₁和θ₂上系统的外力为零，故根据拉格朗日方程即可确定下式
$$\{\begin{array}{l}\frac{d}{dt}\frac{\partial V}{\partial {\dot{\theta }}_1}-\frac{\partial L}{\partial {\theta }_1}=0\\ \frac{d}{dt}\frac{\partial V}{\partial {\dot{\theta }}_2}-\frac{\partial L}{\partial {\theta }_2}=0\end{array}\text{\hspace{1em}(14)}$$
将式(13)带入(14)并化简后即可得到
$$\begin{array}{r}\{\begin{array}{ll}& (\frac{4}{3}{m}_1+4m_2)l_1{\ddot{\theta }}_1+(m_1+2m_2)\ddot{x}cos{\theta }_1+2m_2{l}_{2}cos({\theta }_1-{\theta }_2){\ddot{\theta }}_2\\ & +2m_2{l}_2{\dot{\theta }}_2^{2}sin({\theta }_2-{\theta }_1)+(m_1+2m_2)gsin{\theta }_1=0\\ & \frac{4}{3}{l}_2{\ddot{\theta }}_2-\ddot{x}cos{\theta }_2+2l_1\ddot{x}cos\left({\theta }_1-{\theta }_2\right)+2l_1{\dot{\theta }}_1^{2}sin\left({\theta }_1-{\theta }_2\right)-gsin{\theta }_2=0\end{array}\end{array}\text{\hspace{1em}(15)}$$
为了进一步的化简，作以下设定：
$$\{\begin{array}{l}{\ddot{\theta }}_1=f_1(x,{\theta }_1,{\theta }_2,\dot{x},{\dot{\theta }}_1,{\dot{\theta }}_2)=0\\ {\ddot{\theta }}_2=f_2(x,{\theta }_1,{\theta }_2,\dot{x},{\dot{\theta }}_1,{\dot{\theta }}_2)=0\end{array}\text{\hspace{1em}(16)}$$
之后对各个状态变量，进行一阶展开即可得到公式(17)
$$\{\begin{array}{l}{\ddot{\theta }}_1=k_{11}x+k_{12}{\theta }_1+k_{13}{\theta }_2+k_{14}\dot{x}+k_{15}{\dot{\theta }}_1+k_{16}{\dot{\theta }}_2+k_{17}\ddot{x}\\ {\ddot{\theta }}_2=k_{21}x+k_{22}{\theta }_1+k_{23}{\theta }_2+k_{24}\dot{x}+k_{25}{\dot{\theta }}_1+k_{26}{\dot{\theta }}_2+k_{27}\ddot{x}\end{array}\text{\hspace{1em}(17)}$$
其中
$$\begin{array}{rlr}& k_{11}=\frac{\partial f_1}{\partial x}{|}_{\left(x,{\theta }_1,{\theta }_2,\dot{x},{\theta }_1,{\theta }_2,\ddot{x}\right)=0}=0& k_{21}=\frac{\partial f_2}{\partial x}{|}_{\left(x,{\theta }_1,{\theta }_2,\dot{x},{\theta }_1,{\theta }_2,\ddot{x}\right)=0}=0\\ & k_{12}=\frac{\partial f_1}{\partial {\theta }_1}{|}_{\left(x,{\theta }_1,{\theta }_2,\dot{x},{\theta }_1,{\theta }_2,\ddot{x}\right)=0}=\frac{3m_1+6m_2}{4m_1+3m_2}& k_{22}=\frac{\partial f_2}{\partial {\theta }_1}{|}_{\left(x,{\theta }_1,{\theta }_2,\dot{x},{\theta }_1,{\theta }_2,\ddot{x}\right)=0}=\frac{-9(m_1+2m_2)}{2(4m_1+3m_2)}\frac{g}{l2}\\ & k_{13}=\frac{\partial f_1}{\partial {\theta }_2}{|}_{\left(x,{\theta }_1,{\theta }_2,\dot{x},{\theta }_1,{\theta }_2,\ddot{x}\right)=0}=\frac{-9m_2}{8m_1+6m_2}\frac{g}{l1}& k_{23}=\frac{\partial f_2}{\partial {\theta }_2}{|}_{\left(x,{\theta }_1,{\theta }_2,\dot{x},{\theta }_1,{\theta }_2,\ddot{x}\right)=0}=\frac{3(m_1+3m_2)}{2(2m_1+3m_2)}\frac{g}{l2}\\ & k_{14}=\frac{\partial f_1}{\partial \dot{x}}{|}_{\left(x,{\theta }_1,{\theta }_2,\dot{x},{\theta }_1,{\theta }_2,\ddot{x}\right)=0}=0& k_{24}=\frac{\partial f_2}{\partial \dot{x}}{|}_{\left(x,{\theta }_1,{\theta }_2,\dot{x},{\theta }_1,{\theta }_2,\ddot{x}\right)=0}=0\\ & k_{15}=\frac{\partial f_1}{\partial {\dot{\theta }}_1}{|}_{\left(x,{\theta }_1,{\theta }_2,\dot{x},{\theta }_1,{\theta }_2,\ddot{x}\right)=0}=0& k_{25}=\frac{\partial f_2}{\partial {\dot{\theta }}_1}{|}_{\left(x,{\theta }_1,{\theta }_2,\dot{x},{\theta }_1,{\theta }_2,\ddot{x}\right)=0}=0\\ & k_{16}=\frac{\partial f_1}{\partial {\dot{\theta }}_2}{|}_{\left(x,{\theta }_1,{\theta }_2,\dot{x},{\theta }_1,{\theta }_2,\ddot{x}\right)=0}=0& k_{26}=\frac{\partial f_2}{\partial {\dot{\theta }}_2}{|}_{\left(x,{\theta }_1,{\theta }_2,\dot{x},{\theta }_1,{\theta }_2,\ddot{x}\right)=0}=0\\ & k_{17}=\frac{\partial f_1}{\partial \ddot{x}}{|}_{\left(x,{\theta }_1,{\theta }_2,\dot{x},{\theta }_1,{\theta }_2,\ddot{x}\right)=0}=\frac{3(m_1+m_2)}{4(2m_1+3m_2)}& k_{27}=\frac{\partial f_2}{\partial \ddot{x}}{|}_{\left(x,{\theta }_1,{\theta }_2,\dot{x},{\theta }_1,{\theta }_2,\ddot{x}\right)=0}=\frac{3m_1}{2l_2(4m_1+3m_2)}\end{array}$$

关于这里的线性化方法，大家最好不要用手算，可以去了解以下用matlab的泰勒级数展开来计算，或者对应求偏导，对应函数为 taylar和partial
同时我这里的倒立摆的参数给的是不完整的参数，完整的参数如下，其中m3为转轴重量

$$\begin{array}{rlr}& k_{12}=\frac{\partial f_1}{\partial {\theta }_1}{|}_{\left(x,{\theta }_1,{\theta }_2,\dot{x},{\theta }_1,{\theta }_2,\ddot{x}\right)=0}=\frac{3(m_1+2m_2+2m_3)g}{(4m_1+3m_2+12m_3)l_1}& k_{22}=\frac{\partial f_2}{\partial {\theta }_1}{|}_{\left(x,{\theta }_1,{\theta }_2,\dot{x},{\theta }_1,{\theta }_2,\ddot{x}\right)=0}=\frac{-9(m_1+2m_2+2m_3)g}{2(4m_1+3m_2+12m_3)l_2}\\ & k_{13}=\frac{\partial f_1}{\partial {\theta }_2}{|}_{\left(x,{\theta }_1,{\theta }_2,\dot{x},{\theta }_1,{\theta }_2,\ddot{x}\right)=0}=\frac{-9m_{2}g}{2(4m_1+3m_2+12m_3)l_1}& k_{23}=\frac{\partial f_2}{\partial {\theta }_2}{|}_{\left(x,{\theta }_1,{\theta }_2,\dot{x},{\theta }_1,{\theta }_2,\ddot{x}\right)=0}=\frac{3(m_1+3m_2+3m_3)g}{(4m_1+3m_2+12m_3)l_2}\\ & k_{17}=\frac{\partial f_1}{\partial \ddot{x}}{|}_{\left(x,{\theta }_1,{\theta }_2,\dot{x},{\theta }_1,{\theta }_2,\ddot{x}\right)=0}=\frac{3(2m_1+m_2+4m_3)}{2(4m_1+3m_2+12m_3)l_1}& k_{27}=\frac{\partial f_2}{\partial \ddot{x}}{|}_{\left(x,{\theta }_1,{\theta }_2,\dot{x},{\theta }_1,{\theta }_2,\ddot{x}\right)=0}=\frac{3(m_1-m_3{l}_2)}{2(4m_1+3m_2+12m_3)l_2}\end{array}$$

将x、θ₁、θ₂、x、dθ₁、dθ₂作为系统的状态变量，ddx为控制输入（控制小车的力作为控制输入），则系统的状态空间表达式可写为
$$\{\begin{array}{l}\dot{X}=AX+Bu\\ Y=CX+Du\end{array}\text{\hspace{1em}(18)}$$
其中
$$X=\left[\begin{array}{c}x\\ {\theta }_1\\ {\theta }_2\\ \dot{x}\\ {\dot{\theta }}_1\\ {\dot{\theta }}_2\end{array}\right],Y=\left[\begin{array}{c}x\\ {\theta }_1\\ {\theta }_2\end{array}\right],A=\left[\begin{array}{cccccc}0& 0& 0& 1& 0& 0\\ 0& 0& 0& 0& 1& 0\\ 0& 0& 0& 0& 0& 1\\ 0& 0& 0& 0& 0& 0\\ 0& k_{12}& k_{13}& 0& 0& 0\\ 0& k_{22}& k_{23}& 0& 0& 0\end{array}\right]$$
$$B=\left[\begin{array}{c}0\\ 0\\ 0\\ 1\\ k_{17}\\ k_{27}\end{array}\right],C=\left[\begin{array}{cccccc}1& 0& 0& 0& 0& 0\\ 0& 1& 0& 0& 0& 0\\ 0& 0& 1& 0& 0& 0\end{array}\right],D=\left[\begin{array}{c}0\\ 0\\ 0\end{array}\right]$$
此即系统在平衡点处的线性化模型

线性化模型的第二种方法

由于在平衡点附近运动时，摆动角度很小，所以可以近似认为
$$\begin{gathered} sin\theta =\theta \\ cos\theta =1 \end{gathered}$$
在这种近似条件下，带入之前的拉格朗日方程中，也可以求出对应的线性化模型

详细的求解方式，大家可以参照文献[1]

总结

使用泰勒展开的方式求得线性化模型相比于近似的方法更偏理想化

参考文献
[1] 张立迎. 直线二级倒立摆稳定控制研究[D].山东大学,2009.
[2] 史晓霞,张振东,李俊芳,杨屹.二级倒立摆系统数学模型的建立及意义[J].河北工业大学学报,2001(05):48-51.
[3] 姜辉. 二级倒立摆控制方法研究[D].哈尔滨工程大学,2010.

如果觉得这篇文章对你有帮助，可以在下方为我点个赞，你的支持就是我写作的动力哦
如如果觉得非常好，也可以点击打赏赏博主一点点小钱钱哦

下一篇讲解如何使用simscape对二阶倒立摆对系统进行建模