Matlab 仿真——单自由度倒立摆(1)系统建模

232
0
2021年1月19日 09时03分

文章目录

 

    • 1. 受控对象与设计要求
    • 2. 力分析与系统方程
      • 2.1 转换方程
      • 2.2 状态空间
    • 3. Matlab表达
      • 3.1 转换方程
      • 3.2 状态空间
    • 4. 引用

1. 受控对象与设计要求

 

该例的系统包含一个装有一个倒立摆的小车。我们的控制目标是通过给小车作用一个力,使顶部的倒立摆不落下来。下图标示出了各个变量的含义。

 

在这里插入图片描述

 

对于这个例子,我们有以下参数:

 

(M) 小车质量 0.5 kg

(m) 倒立摆质量 0.2 kg

(b) 小车的摩擦系数 0.1 N/m/sec

(l) 转动点到倒立摆质心的长度 0.3 m

(I) 倒立摆的转动惯量 0.006 kg.m^2

(F) 小车受到的力

(x) 小车的位置坐标

(theta) 倒立摆与垂线的角度

 

我们会用不同的方法设计该系统的控制器:PID,根轨迹,频域分析,状态空间。由于PID,根轨迹以及频域分析法最适合SISO系统的分析研究,因此用这三种方法的时候我们不考虑小车的位置,仅仅考虑如何控制好倒立摆的角度。我们会设计一款控制器,使得小车收到一个冲击(1Nsec)时仍能保持倒立摆垂直,设计要求是5s内稳定倒立摆,且倒立摆角度与稳定状态下的角度相比不超过0.05 弧度。倒立摆的初始位置垂直于地面。因此,我们可以总结出我们的设计需求:

 

  1. θ的稳定时间 < 5s
  2. |θ-θ0| < 0.05 radians

 

而对于状态空间设计方法,我们得以能够处理多输出系统。因此,在用状态空间设计法的时候我们尝试同时控制倒立摆的角度以及小车的位置。以下是我们的设计需求:

 

  1. x 与 θ 的稳定时间 < 5s
  2. x 的上升时间 < 0.5s
  3. |θ-θ0| < 0.05 radians (也就是20°)
  4. 对于x和θ来说,稳态误差 < 2%

2. 力分析与系统方程

 

下图是该系统的受力分析图:

在这里插入图片描述

 

水平方向对小车进行受力分析得到(1)式:

 

捕获

 

你也可以在竖直方向对小车进行受力分析,但没什么卵用。接着在水平方向对倒立摆进行受力分析,得到(2)式 (这一步推导有点快,请参考这个文档):

 

捕获1

 

接着把(2)式代入到(1)式得到该系统的第一个动力学公式:

 

第一个动力学方程

 

捕获2

 

为了得到第二个动力学公式,我们在垂直于倒立摆的方向上列出力平衡方程:

 

捕获3

 

左边的三项不难理解,右边的第一项为引起转动加速度的力,第二项为引起平动加速度的力。为了消去P和N,我们再列出关于倒立摆重心的转动惯量方程:

 

捕获

 

稍微整合一下上述两个方程得到我们要用的第二个动力学方程:

 

第二个动力学方程

 

捕获1

 

上述两个动力学方程是非线性的,为了采用线性系统的方法去设计控制器,还需要把上述两个动力学方程线性化一下。因此我们假设倒立摆的平衡位置 θ  = π ,且倒立摆只在平衡位置下小范围内摆动。我们用ϕ \phiϕ 代表倒立摆与平衡位置的角度偏差,那么倒立摆在任意时刻的位置 θ π  + ϕ 。有了上述线性化假设之后,我们可以得到:

 

  1. 捕获2
  2. 捕获3
  3. 捕获4

 

将上述结果带入到我们的动力学方程组中得到:

 

第一个动力学方程

 

捕获5

 

第二个动力学方程

 

捕获6

2.1 转换方程

 

对动力学方程组分别做拉氏变换

 

捕获7

 

要得到输入U ( s ) U(s)U(s)与输出 Φ ( s ) \Phi(s)Φ(s) 的关系,我们需要消掉X ( s ) X(s)X(s),于是经过一系列的化简代入我们最终得到:

 

捕获

 

其中:

 

捕获1

2.2 状态空间

 

有了线性化之后的动力学方程组,一顿化简得到:

 

捕获

 

该系统的输出有两个,一是小车的位置x xx ,二是倒立摆的角度θ \thetaθ ,因此y yy有两项

3. Matlab表达

 

现在我们就可以在Matlab里面表达出我们的系统了

 

3.1 转换方程

 

M = 0.5;
m = 0.2;
b = 0.1;
I = 0.006;
g = 9.8;
l = 0.3;
q = (M+m)*(I+m*l^2)-(m*l)^2;
s = tf('s');

P_cart = (((I+m*l^2)/q)*s^2 - (m*g*l/q))/(s^4 + (b*(I + m*l^2))*s^3/q - ((M + m)*m*g*l)*s^2/q - b*m*g*l*s/q);

P_pend = (m*l*s/q)/(s^3 + (b*(I + m*l^2))*s^2/q - ((M + m)*m*g*l)*s/q - b*m*g*l/q);

sys_tf = [P_cart ; P_pend];

inputs = {'u'};
outputs = {'x'; 'phi'};

set(sys_tf,'InputName',inputs)
set(sys_tf,'OutputName',outputs)

sys_tf

 

输出:

 

在这里插入图片描述

3.2 状态空间

 

M = .5;
m = 0.2;
b = 0.1;
I = 0.006;
g = 9.8;
l = 0.3;

p = I*(M+m)+M*m*l^2; %denominator for the A and B matrices

A = [0      1              0           0;
     0 -(I+m*l^2)*b/p  (m^2*g*l^2)/p   0;
     0      0              0           1;
     0 -(m*l*b)/p       m*g*l*(M+m)/p  0];
B = [     0;
     (I+m*l^2)/p;
          0;
        m*l/p];
C = [1 0 0 0;
     0 0 1 0];
D = [0;
     0];

states = {'x' 'x_dot' 'phi' 'phi_dot'};
inputs = {'u'};
outputs = {'x'; 'phi'};

sys_ss = ss(A,B,C,D,'statename',states,'inputname',inputs,'outputname',outputs)

 

输出

 

在这里插入图片描述

我们可以通过把状态空间表达转换为转换方程:

 

sys_tf = tf(sys_ss)

 

输出

 

在这里插入图片描述

 

我们发现上述转换方程里面有系数非常接近于零,这是由于Matlab在四舍五入的时候带来的小误差。我们完全可以手动把这些系数改为0。完了之后我们这里转换得到的转换方程和之前我们直接得到转换方程表征的应该是完全一样的系统

 

4. 引用

 

https://ctms.engin.umich.edu/CTMS/index.php?example=InvertedPendulum&section=SystemModeling

发表评论

后才能评论