一级倒立摆系统分析.doc

一级倒立摆的系统分析
一、 倒立摆系统的模型建立
如图1-1所示为一级倒立摆的物理模型

小 车
导 轨
x
l
F
摆杆
ɸɸ

图1-1 一级倒立摆物理模型
对于上图的物理模型我们做以下假设：
M：小车质量
m：摆杆质量
b：小车摩擦系数
l：摆杆转动轴心到杆质心的长度
I：摆杆惯量
F：加在小车上的力
x：小车位置
ɸ：摆杆与垂直向上方向的夹角
θ：摆杆与垂直向下方向的夹角（考虑到摆杆初始位置为竖直向下）
图1-2是系统中小车和摆杆的受力分析图。其中，N和P为小车与摆杆相互作用力的水平和垂直方向的分量。注意：实际倒立摆系统中的检测和执行装置的正负方向已经完全确定，因而矢量方向定义如图所示，图示方向为矢量正方向。
F
N
P
M
X
Bx’
X’’
mg
N
P
θ
Iθ’’

图1-2 小车及摆杆受力分析
分析小车水平方向受力，可以得到以下方程：
Mx＝Ｆ-ｂx-Ｎ (1-1)
由摆杆水平方向的受力进行分析可以得到以下方程：
N=md2dt2(x+lsinθ) (1-2)
即： N=mx+mlθcosθ-mlθ2sinθ (1-3)
将这个等式代入式（1-1）中，可以得到系统的第一个运动方程：
M+mx+bx+mlθcosθ-mlθ2sinθ=F (1-4)
为推出系统的第二个运动方程，我们对摆杆垂直方向上的合力进行分析，可以得出以下方程：
P-mg=md2dt2(lcosθ) (1-5)
P-mg=- mlθsinθ-mlθ2cosθ (1-6)
利用力矩平衡方程可以有：
-Plsinθ-Nlcosθ=Iθ (1-7)
注意：此方程中的力矩方向，由于θ=π+ɸ，cosɸ=-cosθ，sinɸ=-sinθ，所以等式前面含有负号。
合并两个方程，约去P和N可以得到第二个运动方程：
I+ml2θ+mglsinθ=-mlxcosθ (1-8)
设θ=π+ɸ，假设ɸ与1（单位是弧度）相比很小，即ɸ＜＜1，则可以进行近似处理：cosθ=-1，sinθ=-ɸ，（dθdt）2=0。用u来代表被控对象的输入力F，线性化后的两个运动方程如下：
I+ml2ɸ-mglɸ=mlxM+mx+bx-mlɸ=u (1-9)
假设初始条件为0，则对式（1-9）进行拉普拉斯变换，可以得到：
I+ml2Φss2-mglΦs=mlX(s)s2M+mXss2+bXss-mlΦss2=U(s) (1-10)
由于输出为角度ɸ，求解方程组的第一个方程，可以得到：
Xs=[I+ml2ml-gs2]Φs (1-11)
或改写为：ΦsXs=mls2I+ml2s2-mgl (1-12)
如果令v=x,则有：ΦsV（s）=mlI+ml2s2-mgl (1-13)
如果将上式代入方程组的第二个方程，可以得到：
M+mI+ml2ml-gsΦss2+bI+ml2ml+gs2Φss-mlΦss2=U(s) (1-14)
整理后可得传递函数：
ΦsU(s)=mlqs2s4+bI+ml2qs3-M+mmglqs2-bmglqs (1-15)
其中 q=[M+mI+ml2-(ml)2]
假设系统状态空间方程为：
X=AX+Bu
y=CX+Du (1-16)
方程组对x,ɸ解代数方程，可以得到解如下：
x=xx=-I+ml2bIM+m+Mml2x+m2gl2IM+m+Mml2ɸ+I+ml2IM+m+Mml2uɸ=ɸɸ=-mlbIM+m+Mml2x+mgl(M+m)IM+m+Mml2ɸ+mlIM+m+Mml2u (1-17)
整理后可以得到系统状态空间方程：
xxɸɸ=01000-I+ml2bIM+m+Mml2m2gl2IM+m+Mml2000010-mlbIM+m+Mml2mgl(M+m)IM+m+Mml20xxɸɸ+0I+ml2IM+m+Mml20mlIM+m+Mml2u
y=xɸ=10000010xxɸɸ+00u (1-18)
由（1-9）的第一个方程为：
I+ml2ɸ-mgl ɸ=mlx
对于质量均匀分布的摆杆可以有：
I=13ml2
于是可以得到：13ml2+ml2ɸ-mgl ɸ=mlx
化简可以得到：ɸ=3g4lɸ+34lx (1-19)
设X={x, x, ɸ , ɸ}，u=x则有：
xx ɸɸ=010000000001003g4l0xxɸɸ+01034lu
y=xɸ=10000010xxɸɸ+00u (1-20)
以上公式推理是根据牛顿力学的微分方程验证的。在实际系统中模型参数如下：
M 小车质量 1.096 Kg
m 摆杆质量 0.109 Kg
b 小车摩擦系数 0 .1N/m/sec
l 摆杆转动轴心到杆质心的长度 0.2 5m
I 摆杆惯量 0.0034 kg*m*m
将上述参数代入，就可以得到系统的实际模型。
摆杆角度和小车位移的传递函数：
ΦsXs=0.02725s20.0102125s2-0.26705 (1-21)
摆杆角度和小车加速度之间的传递函数为：
ΦsV（s）=0.027250.0102125s2-0.26705(1-22)
摆杆角度和小车所受外界作用力的传递函数：
ΦsU（s）=2.35655ss3+0.0883167s2-27.9169s-2.30942(1-23)
以外界作用力作为输入的系统状态方程：
xx ɸɸ=01000-0.08831670.629317000010-0.23565527.82850xxɸɸ+00.88316702.35655u
y=xɸ=10000010xxɸɸ+00u (1-24)
以小车加速度作为输入的系统状态方程：
xx ɸɸ=0100000000010029.40xxɸɸ+0103u
y=xɸ=10000010xxɸɸ+00u （1-25）
综述可知以上就是一级倒立摆系统的模型建立过程，最终得出了实际模型的传递函数和状态空间方程。
二、 系统模型的转换
以小车加速度作为输入的系统状态方程为例，将系统状态方程转化为能控标准型，能观标准型和约当标准型。
由系统状态方程可知：
A=0100000000010029.40
B=0103
C=10000010
D=00
1、 转化为能控标准型
定出系统特征多项式：
>> a=poly(A) a =
1.0000 -0.0000 -29.4000 0 0
由此可知a0=0, a1=0, a2=-29.4, a3=0。
>> b3=C*B
b3 =
0
0
b2=C*A*B+a3*C*B
b2 =
1
3
>> b1=C*A^2*B+a3*C*A*B+a2*C*B
b1 =
0
0
>> b0=C*A^3*B+a3*C*A^2*B+a2*C*A*B+a1*C*B
b0 =
-29.4000
0
所以系统的能控标准型为：
x1x1 ɸ1ɸ1=0 10000 1000010029.40x1x1ɸ1ɸ1+ 0001u
y=-29.40100030x1x1ɸ1ɸ1+00u
2、 转化为能观标准型
利用对偶性求出能观标准型为：
x1x1 ɸ1ɸ1=0 00010 0001029.40010x1x1ɸ1ɸ1+ 0-29.4003300u
y=[0001]x1x1ɸ1ɸ1
3、 转化为约当标准型
首先求出系统的特征值以及相应的特征向量：
A=[0 1 0 0;0 0 0 0;0 0 0 1;0 0 29.4 0]
A =
0 1.0000 0 0
0 0 0 0
0 0 0 1.0000
0 0 29.4000 0
>> [V,D]=eig(A)
V =
0 0 1.0000 -1.0000
0 0 0 0.0000
0.1814 -0.1814 0 0
0.9834 0.9834 0 0
D =
5.4222 0 0 0
0 -5.4222 0 0
0 0 0 0
0 0 0 0
其中D表示A全部特征值构成的对角阵，V表示相对应的特征向量。
求出变换矩阵V的逆：
>> V1=inv(V)
Warning: Matrix is close to singular or badly scaled.
Results may be inaccurate. RCOND = 1.720635e-292.
V1 =
1.0e+291 *
0 0 0.0000 0.0000
0 0 -0.0000 0.0000
0.0000 2.4948 0 0
0 2.4948 0 0
计算变换后的系数矩阵：
>> A1=V1*A*V
A1 =
5.4222 0 0 0
0.0000 -5.4222 0 0
0 0 0 0.0000
0 0 0 0
B1=V1*B

B1 =
1.0e+291 *
0.0000
0.0000
2.4948
2.4948
所以系统的约当标准型为：
x1x1 ɸ1ɸ1=5.4222 0000-5.4222 0000000000x1x1ɸ1ɸ1+ 1.0e+291 *002.49482.4948u
三、 开环阶跃响应曲线及分析
利用已知的状态空间方程来进行阶跃响应分析，在MATLAB中可以写入以下命令：
>> A=[0 1 0 0;0 0 0 0;0 0 0 1;0 0 29.4 0]；
>> B=[0;1;0;3]；
>> C=[1 0 0 0;0 1 0 0]；
>> D=[0;0]；
>> step(A,B,C,D)

可以看出，在单位阶跃响应作用下，小车位置和摆杆角度都是发散的。
四、 判断系统稳定性
判断系统的稳定性可以利用根轨迹来判断，已知实际系统的开环传递函数为：ΦsV（s）=0.027250.0102125s2-0.26705，则其根轨迹图形可以利用MATLAB键入如下命令来完成。
>> num=[0.02725];
>> den=[0.0102125 0 -0.26705];
>> z=roots(num)
z =
Empty matrix: 0-by-1
>> p=roots(den)
p =
5.1136
-5.1136
>> rlocus(num,den)

可以看出系统没有零点，有两个极点，并且有一个极点为正。由画出的根轨迹图形可以看出闭环传递函数的一个极点位于复平面的右半平面，这就意味着系统是一个不稳定的系统。
五、 能控性和能观性分析
对于系统的能控性和能观性分析，可以利用能控性秩判据和能观性秩判据。
能控性秩判据：对于n维连续时间线性时不变系统，构成能控性判别矩阵：Qc=[B⋮AB⋮⋯⋮An-1B]，则系统完全能控的充要条件为：
rankQc=rankB⋮AB⋮⋯⋮An-1B=n
能观性秩判据：对于n维连续时间线性时不变系统，构成能观性判别矩阵：Qo=CCA⋮CAn-1，则系统完全能观的充要条件为：
rankQo=rankCCA⋮CAn-1=n
利用MATLAB键入以下命令来进行判断：
>> A=[0 1 0 0;0 0 0 0;0 0 0 1;0 0 29.4 0];
>> B=[0;1;0;3];
>> C=[1 0 0 0;0 1 0 0];
>> D=[0;0];
>> Qc=[B A*B A^2*B A^3*B]
Qc =
0 1.0000 0 0
1.0000 0 0 0
0 3.0000 0 88.2000
3.0000 0 88.2000 0
>> R1=rank(Qc)
R1 =
4
>> Qo=[C;C*A;C*A^2;C*A^3]
Qo =
1 0 0 0
0 1 0 0
0 1 0 0
0 0 0 0
0 0 0 0
0 0 0 0
0 0 0 0
0 0 0 0
>> R2=rank(Qo)
R2 =
2
可以看出，系统的完全能控矩阵的秩等于系统的状态变量维数，系统的输出完全能观测矩阵的秩等于系统输出向量y的维数，所以系统是可以完全能控完全能观测的系统。
六、 根轨迹校正以及仿真
已知系统的传递函数：
Gs=0.027250.0102125s2-0.26705
设计控制器使得调整时间ts=0.5s(2%)；最大超调Mp≤10%。
计算整理可得超前校正装置的零点和极点分别为：
zc=-6.92214；zp=-26.4568，由此可得校正后的传递函数：
Q=GsKs=K(s+6.92214)s+26.45680.027250.0102125s2-0.26705
利用MATLAB命令观察校正后的根轨迹图形：

图中可以看出，系统三条根轨迹都位于左半平面，只需要选定适当的K就可以使得系统稳定。
现在利用MATLAN仿真，选择不同的K观察系统的阶跃响应曲线。
K=141.137时，画出闭环控制系统：

响应曲线图：

K=120.115时，响应曲线：

K=50.129时的响应曲线：

K=10.110时的响应曲线：

K=1.137时的响应曲线：

对比不同的K值变化可知，闭环控制系统是稳定的。