旋转倒立摆系统设计分析报告.docx

旋转倒立摆控制系统建模、分析与设计
李振
(1. 上海大学机电工程与自动化学院 上海 200072)

摘要：本文介绍了单级旋转倒立摆系统的构成，并利用Lagrange方程建立了数学模型，根据倒立摆系统的数学模型在Matlab环境下设计了控制器，并进行仿真。设计了LQR状态跟随运动控制器。结果表明，LQR状态跟随器控制系统超调小，响应速度快，鲁棒性强。根据设定的Q和R，计算出LQR最优控制参数K。
关键词：旋转倒立摆；状态反馈；最优控制策略LQR
中图分类号：TG214

Modeling Analysis and Design of Rotary Inverted Pendulum Control System
LI Zhen
(1.College of Mechatronic Engineering and automation, Shanghai University, Shanghai 200072)

Abstract：The structure andmathematic model of Single-Rotational Inverted Pendulum(RIP) is presented in this paper. The mathematical model is set on the Language equation. According to the mathematical model of inverted pendulum system controller was designed in Matlab environment, and carry on the simulation. LQR state following the motion controller is designed. The results show that the LQR state follower control system has small overshoot, fast response and strong robustness. According to the set of Q and R, the LQR optimal control parameter K is calculated.
Key words：Rotational Inverted Pendulum State-feedback Controller LQR controller
0 前言
常见的倒立摆系统主要包括直线倒立摆、旋转型倒立摆和环型倒立摆。他们之间主要是机械结构的不同，原理上大同小异。各种倒立摆系统在建模方法、模型分析上十分类似，而其主要区别在于采取的控制理论不同。
旋转倒立摆的摆杆旋臂是在竖直平面内旋转的，一般的结构如图1所示。

图1 旋转倒立摆的总体结构
在本实验中，为了建模过程方便，对旋转倒立摆的结构进行了简化。利用现代控制理论的有关知识，对简化的旋转倒立摆系统建立状态空间模型，分析其能控及能观性并进行状态反馈控制的设计。
1 旋转倒立摆控制系统建模
本试验仅对旋转倒立摆的简易模型进行设计与分析。
假设旋臂仅在水平面内运动，忽略空气阻力及摩擦，得到简易旋转倒立摆的简化模型如图2所示。该模型描述了旋转倒立摆的运动过程。由于旋臂在电机的带动下旋转，摆杆随之摆动，通过对旋臂的旋转控制，保持摆杆在竖直平面内。

图2 旋转倒立摆的简化模型
建模过程中，机械结构参数有：①摆杆的集中质量m②旋臂长度r③摆杆等效长度L。
设摆杆偏离垂直向上的角度为α，旋臂转过的角度为θ。
该系统的状态空间方程推导如下：
将摆杆竖直面的运动在x与y方向上分解：
Vm=Lcosαα+Lsinαα (1)
旋臂的速度：
VM=rθ (2)
则有摆杆实际运动的速度分解为：
Vx=rθ-Lcosαα (3)
Vy=-Lsinαα (4)
系统的总势能为重力势能：

EP=mgh=mgLcosα (5)
系统总动能的表达式：
Ek=EkM+EkVx+EkVy+Ekm (6)
求摆杆的转动惯量：
I=x2dm=mL-hL-hx2dx=13mL2-3Lh+3h2
代入h=L，得
Jm=13mL2 (7)
则总动能为：
Ek=12JMθ2+12mrθ-Lcosαα2+12m-Lsinαα2+12Jmα2 (8)
定义Lagrange算子La：
La=Ek-Ep
=12JMθ2+23mL2α2-mLrcosααθ+12mr2θ2-mgLcosα (9)
由Lagrange方程：
∂∂t∂La∂θ-∂La∂θ=T-Beqθ (10)
∂∂t∂La∂α-∂La∂α=0 (11)
得：
JMθ-mLrcosαα+mLrsinαα2+mrθ=T-Beqθ (12)
43mL2α-mLrθcosα-mgLsinα=0 (13)
近似认为α很小，对方程(10)、(11)进行线性化，即认为 sinα→α,cosα→1 ，高次项α2→0，得：
JM+mr2θ-mLrα=T-0 (14)
43mL2α-mLrθ-mgLα=0 (15)
其中电机输出力矩为：
T= (16)
其中ηm、ηg为电机、减速器效率；Ki、Kg为电机力矩和减速比系数；Vm为电机电枢电压；Rm为电机电枢电阻；Km为电机反电势系数。
整理上面的式子，可得如下状态空间：
x=Ax+Bu
Y=Cx+Du
其中
x=θθαα Y=θα u=Vm
A=0100-4ηmηgKiKg2KmRm4JM+mr2-4Beq4JM+mr23mgr4JM+mr200001-3rηmηgKiKg2KmLRm4JM+mr2-3rBeqL4JM+mr2gJM+mr2L4JM+mr20
B=04ηmηgKiKgRm4JM+mr203rηmηgKiKgLRm4JM+mr2 C=10000010 D=00
即
θθαα=0100-4ηmηgKiKg2KmRm4JM+mr2-4Beq4JM+mr23mgr4JM+mr200001-3rηmηgKiKg2KmLRm4JM+mr2-3rBeqL4JM+mr2gJM+mr2L4JM+mr20θθαα+04ηmηgKiKgRm4JM+mr203rηmηgKiKgLRm4JM+mr2Vm (17)
θα=10000010θθαα (18)
系统建模完成。
2 系统的性能分析
方程中的一些参数如下表：
表1 方程中的参数
符号
物理意义/单位
数值
L
摆杆质心到转轴的距离/m
0.165
m
摆杆质量/kg
0.105
r
旋臂长度/m
0.254
JM
旋臂的转动惯量/kg·m2
2.0×10-3
g
重力加速度/m·s-2
9.8
Ki
电机力矩系数/N·m·A-1
7.767×10-3
Kg
变速器齿轮比
5:1

Km
反向电势系数/V·s·rad-1
7.767×10-3
Rm
直流电机电枢电阻/Ω
2.6
ηm
直流电机效率/%
69
ηg
变速器效率/%
90
Beq
粘性阻尼系数N·s·m -1
4.0×10-3
将以上数据代入状态空间，得
θθαα=0100-0.0975-1.083053.072200001-0.1126-1.250335.27320θθαα+02.511302.8994Vm (19)
2.1系统能控性分析
a) 能控性定义
线性定常系统的状态方程为
x=Ax+Bu
给定系统一个初始状态x(t0)，如果在t1>t0的有限时间区间[t0, t1]内，存在容许控制u(t)，使x(t1)=0，则称系统状态在t0时刻是能控的；如果系统对任意一个初始状态都能控，则称系统是状态完全能控的。
b) 能控性判据
线性定常系统为状态能控的充分必要条件是以下的矩阵满秩：
QC=[BABA2BA3B…An-1B]
Rank QC=n
轶判据本身很简单，因此是判断能控性最为常用的方法。
c) 系统能控性判断
由系统能控性判据（轶判据），计算能控性矩阵为
Qc=BABA2BA3B=02.5113-2.7197156.57832.5113-2.7197156.5783-335.949602.8994-3.1400105.38922.8994-3.1400105.3892-306.2266 (20)
能控性矩阵的秩 rank Qc = rank A = 4，因此系统能控。
2.1系统能观性分析
a) 能观性定义
线性定常系统方程为x=Ax+Buy=Cx
如果在有限时间区间[t0, t1](t1>t0)内，通过观测y(t)，能够惟一地确定系统的初始状态x(t0)，称系统状态在t0是能观测的。如果对任意的初始状态都能观测，则称系统是状态完全能观测的。
b) 能观性判据
线性定长系统为能观测的充分必要条件是以下能观性矩阵满秩：
rank Q0=n
Q0=CCA⋮CAn-1nm×n
轶判据本身很简单，因此是判断能观性最为常用的方法。
c) 系统能观性判断
由系统能观性判据（轶判据），计算能控观矩阵为
Q0=CCACA2CA3=1000001001000001-0.0975-1.083053.07220-0.1126-1.250335.273200.10561.0753-57.475653.07220.12191.2415-66.385235.2732 (21)
能观性矩阵的秩 rank Q0 = rank A = 4，因此系统能观。
2.1系统稳定性分析
a) 李雅普诺夫判据第一法
线性定常系统方程为x=Ax+Buy=Cx的系统矩阵A的特征根：
A) 全部为负实部根，则实际系统渐进稳定，在线性化过程中忽略的高阶导数项对系统稳定性无影响。即线性化模型渐进稳定，线性化前的非线性化系统也渐进稳定。
B) 只要有一个为正实部根，则实际系统不稳定，在线性化过程中忽略的高阶导数项对系统稳定性无影响。即线性化模型不稳定，线性化前的非线性系统也不稳定。
C) 只要有一个实部为0（纯虚根），其余为负实部，则在线性化过程中忽略的高阶导数项对系统稳定性有影响，不能用线性化模型来判断实际系统的稳定性，必须采用实际系统的非线性模型来判断。
b) 系统稳定性判断
由λ·I-A=0求得系统的特征根为：
λ1= 4.9007, λ2=0.9274 , λ3=-0.0817 , λ4= -6.8294
其中特征根λ1, λ2具有正实部，因此系统不稳定。
3 状态反馈矩阵设计
采用极点配置法设计状态反馈矩阵K，对系统进行状态反馈设计。其中期望的闭环极点可以根据最优线性二次型方法得到。
3.1线性二次型（LQR）方法
最优线性二次型控制的原理为，根据现代控制理论，在系统的状态空间方程
x=Ax+Bu
中，设状态反馈控制律的形式为：ut=K×x(t)，其中K为最佳控制向量矩阵，使线性二次型最优性能指标J达到最小值。
(22)
旋转倒立摆的控制中，LQR算法的控制任务是使价值函数J最小。
J=0∞xTQCx+uTRCudt
(23)
由于这里只有一个变量矩阵，因此可以把上面的函数简化为：
(24)
其中u不受限制，Q是半正定常数矩阵，R为正定常数矩阵。求得状态反馈矩阵K=R-1BTP，P为满足代数Riccati方程的唯一正定对称解。代数方程为：
PA+ATP+Q-PBR-1BTP=0
(25)
LQR控制的原理图为

代数方程是考虑到控制能量的损耗而引进的，加权矩阵Q和R确定了误差和能量损耗的相对重要性。并且假设控制向量u(t)是无约束的。
3.2状态反馈矩阵设计
这里加权矩阵分别为：
Q=2.30004500000000000 R=1
根据期望性能指标选取合适的加权矩阵RC和QC利用MATLAB命令lqr就可以得到最优增益矩阵矩阵K的值。该LQR最优控制可以使系统从任意初始状态返回到平衡状态（0状态）。
K=lqr( A,B,Q,R)
获取的K矩阵为：
K=1.47821.324533.35283.3813
4 结果分析
开环及状态反馈的单位阶跃响应的响应曲线如下图：

图1 系统开环单位阶跃响应（不稳定）

图2 状态反馈系统的单位阶跃响应（稳定）
附录：Matlab仿真程序
clc;
clear all;
close all;
%参数值
L = 0.165; m = 0.105; r=0.254; J_M = 2.0 * 10^(-3);
g = 9.8; K_i = 7.767 * 10^(-3); K_m = 7.767 * 10^(-3);
K_g = 5; R_m = 2.6; ita_m = 0.69; ita_g = 0.90; B_eq = 4.0 * 10^(-3);
%建立状态空间方程
e1= ita_m * ita_g * K_i * K_g;
e2 = e1 * K_g * K_m;
e3 = 4 * J_M + m * r^2 ;
A21 = - (4 * e2 ) / ( R_m * e3 );
A22 = - (4 * B_eq) / e3;
A23 = (3 * m * g * r) / e3;
A41 = - (3 * r * e2 ) / ( L * R_m * e3) ;
A42 = - (3 * r * B_eq) / ( L * e3) ;
A43 = ( (J_M + m * r^2 ) * g) / (L * e3);
B2 = (4 * e1 ) / ( R_m * e3 ) ;
B4 = (3 * r * e1 ) / ( L * R_m * e3) ;
A = [0, 1, 0, 0 ; A21, A22, A23 0 ; 0, 0, 0, 1 ; A41, A42, A43, 0];
B = [0; B2; 0; B4];
C = [1, 0, 0, 0 ; 0, 0, 1, 0];
D = [0; 0];
G=ss(A, B, C, D);
%能观能控性分析
n=size(A,1);
Q_c=B;
for i=1: 1: n-1
H = (A^i) * B;
Q_c=[Q_c , H];
end
if ( rank(Q_c) == n )
display('系统能控')
else
display('系统不能控')
end
Q0=C;
for i=1: 1: n-1
L = C * (A^i) ;
Q0=[Q0 ; L];
end
if ( rank(Q0) == n )
display('系统能观')
else
display('系统不能观')
end
%稳定性分析
[num,den]=ss2tf(A ,B, C, D);
s = roots(den);
k = size(s,1);
for l = 1:1:k
if ( real (s(l) ) > 0 )
display('系统不稳定')
break;
else
if ( l == k)
display('系统稳定')
end
end
end
figure(1)
step(G)
%LQR控制器设计
Q = [2.3, 0, 0 , 0; 45, 0, 0, 0; 0, 0, 0 , 0;0, 0, 0, 0];
R=1;
K=lqr(A, B, Q,R);

%状态反馈设计
a = A-B * K;
G0 = ss(a , B, C, D);
figure(2)
step(G0)
参 考 文 献
[1] Keyu Lai, Jin Xiao, Xiaoguang Hu, Jianxin Fan and Bing Wu, "Modeling and control for stability and rotation velocity of a rotary inverted pendulum," Industrial Electronics and Applications (ICIEA), 2015 IEEE 10th Conference on, Auckland, 2015, pp. 955-960.
[4] M. Akhtaruzzaman and A. A. Shafie, "Modeling and control of a rotary inverted pendulum using various methods, comparative assessment and result analysis," Mechatronics and Automation (ICMA), 2010 International Conference on, Xi'an, 2010, pp. 1342-1347
[5] S. M. M. Anvar, I. Hassanzadeh and G. Alizadeh, "Design and implementation of sliding mode-state feedback control for stabilization of Rotary Inverted Pendulum," Control Automation and Systems (ICCAS), 2010 International Conference on, Gyeonggi-do, 2010, pp. 1952-1957.
[2] 刘二林, 姜香菊. 单级旋转倒立摆运动控制研究[J]. 制造业自动化, 2015(13):48-51.
[3] 段旭东, 许可. 单级旋转倒立摆的建模与控制仿真[J]. 机器人技术与应用, 2002(5):43-46.
[4] 吴爱国, 张小明, 张钊. 基于Lagrange方程建模的单级旋转倒立摆控制[J]. 中国工程科学, 2005(10):11-15.
[5] 张姝, 朱善安. 基于Lagrange方程建模的旋转倒立摆的分析与控制[J]. 工业控制计算机, 2003, 16(8):30-32.
[6] 赵平平. 用MATLAB进行现代控制理论分析与设计[J]. 科技信息:科学教研, 2007(35).