一阶倒立摆控制新版专业系统设计.doc

课程设计说明书

课程名称： 控制系统课程设计
设计题目：一阶倒立摆控制器设计
院 系： 信息和电气工程学院
班 级：
设 计 者：
学 号：
指导老师：
设计时间：2月25日到3月8号

课程设计（论文）任务书

专 业
自动化
班 级
0902101
学 生

指导老师

题 目
一阶倒立摆课程设计
子 题

设计时间
2 月 25 日 至 年 3 月 8 日 共 2 周
设计要求
设计（论文）任务和基础要求，包含设计任务、查阅文件、方案设计、说明书（计算、图纸、撰写内容及规范等）、工作量等内容。

1．建立一阶倒立摆数学模型
2．做模型仿真试验
（1）给出Matlab仿真程序。
（2）给出仿真结果和响应曲线。
3．倒立摆系统PID控制算法设计
设计PID控制器，使适当在小车上施加1N脉冲信号时，闭环系统响应指标为：
（1）稳定时间小于5秒
（2）稳态时摆杆和垂直方向夹角改变小于0.1 弧度
并作PID控制算法MATLAB仿真
4．倒立摆系统最优控制算法设计
用状态空间法设计控制器，使适当在小车上施加0.2m阶跃信号时，闭环系统响应指标为：
（1）摆杆角度和小车位移稳定时间小于5秒
（2）上升时间小于1秒
（3）超调量小于20度（0.35弧度）
（4）稳态误差小于2%。

指导老师签字： 系（教研室）主任签字：
3月 5日
目录

一、 建立一阶倒立摆数学模型 4
1. 一阶倒立摆微分方程模型 4
2. 一阶倒立摆传输函数模型 6
3. 一阶倒立摆状态空间模型 7
二、 一阶倒立摆matlab仿真 9
三、 倒立摆系统PID控制算法设计 13
四、倒立摆系统最优控制算法设计 23
五、 总结 28
六、 参考文件 29

一、 建立一阶倒立摆数学模型
　　　首先建立一阶倒立摆物理模型。在忽略空气阻力和多种摩擦以后， 可将直线一级倒立摆系统抽象成小车和匀质杆组成系统，图1所表示。

　　　系统内部各相关参数定义以下：
　　　M 小车质量
　　　m 摆杆质量
　　　b 小车摩擦系数
　　　l 摆杆转动轴心到杆质心长度
　　　I 摆杆惯量
　　　F 加在小车上力
　　　x 小车位置
　　　φ 摆杆和垂直向上方向夹角
θ 摆杆和垂直向下方向夹角（考虑到摆杆初始位置为竖直向下）
1. 一阶倒立摆微分方程模型
对一阶倒立摆系统中小车和摆杆进行受力分析，其中，N和 P为小车和摆杆相互作用力水平和垂直方向分量。

图 1-2 小车及摆杆受力图

分析小车水平方向所受协力，能够得到以下方程：
（1-1）
由摆杆水平方向受力进行分析能够得到下面等式：
（1-2）
即：
（1-3）

把这个等式代入式(1-1)中，就得到系统第一个运动方程：

（1-4）

为了推出系统第二个运动方程，我们对摆杆垂直方向上协力进行分析，能够得到下面方程：
（1-5）
即：
（1-6）

力矩平衡方程以下：
（1-7）
因为所以等式前面有负号。
合并这两个方程，约去 P和 N，得到第二个运动方程：
（1-8）

设 ，（φ是摆杆和垂直向上方向之间夹角），假设φ <<1弧度， 则能够进行近似处理：。用u代表被控对象输入力F，利用上述近似进行线性化得直线一阶倒立摆微分方程为：
（1-9）

2. 一阶倒立摆传输函数模型
对式(1-9)进行拉普拉斯变换，得：
（2-1）

注意：推导传输函数时假设初始条件为 0。

因为输出为角度φ，求解方程组第一个方程，可得：
（2-2）

或
（2-3）

假如令，则有：
（2-4）

（2-5）
把上式代入方程组（2-1）第二个方程，得：

整理后得到传输函数：
（2-6）

其中。
3. 一阶倒立摆状态空间模型
设系统状态空间方程为：
（3-1）

方程组（2-9）对解代数方程，得到解以下：
（3-1）

整理后得到系统状态空间方程：
（3-2）

（3-3）

摆杆惯量为，代入（1-9）第一个方程为：

得：

化简得：

（3-4）

设 则有：

（3-5）

4.实际系统传输函数和状态方程
实际系统模型参数以下:
M 小车质量 0.5 Kg
m 摆杆质量 0.2 Kg
b 小车摩擦系数 0 .1N/m/sec
l 摆杆转动轴心到杆质心长度 0.3m
I 摆杆惯量 0.006 kg*m*m

代入上述参数可得系统实际模型。
摆杆角度和小车位移传输函数：
（4-1）

（4-2）

摆杆角度和小车加速度之间传输函数为：

（4-3）

摆杆角度和小车所受外界作用力传输函数：

以外界作用力作为输入系统状态方程：

（4-4）

以小车加速度为输入系统状态方程：
（4-5）

二、 一阶倒立摆matlab仿真

实际系统参数以下，根据上面给出例子求系统传输函数、状态空间方程，并进行脉冲响应和阶跃响应matlab仿真。
小车质量 1.096Kg
摆杆质量 0.109Kg
小车摩擦系数 0.1 N/m/sec
摆杆转动轴心到杆质心长度 0.25m
摆杆惯量 0.0034 kg*m*m
采样时间 0.005秒
1. 传输函数法
Matlab程序以下：
M=1.096;
m=0.109;
b=0.1;
I=0.0034;
g=9.8;
L=0.25;
q=(M+m)*(I+m*L^2)-(m*L)^2;
num=[m*L/q 0 0]
den=[1 b*(I+m*L^2)/q -(M+m)*m*g*L/q -b*m*g*L/q 0];
[r,p,k]=residue(num,den);
s=p;
得到传输函数分子：
num =

2.3566 0 0
和传输函数分母:
den =

1.0000 0.0883 -27.8285 -2.3094 0
开环极点：
s =

-5.2780
5.2727
-0.0830
0

由此可知，系统传输函数多项式表示式为：
（2-1）
系统开环极点为（s):、、、，因为有一个开环极点在平面右半部，开环系统并不是稳定。
系统脉冲响应以下，由图也可见，系统并不稳定。

图2.1 开环系统脉冲响应

2. 状态空间法
状态空间法能够进行单输入多输出系统设计，所以在这个试验中，我们将尝试同时对摆杆角度和小车位置进行控制。为了更具挑战性，给小车加一个阶跃输入信号。
我们用 Matlab 求出系统状态空间方程各矩阵，并仿真系统开环阶跃响应。在这里给出一个state.m文件，实施这个文件，Matlab将会给出系统状态空间方程A，B，C和D矩阵，并绘出在给定输入为一个0.2 N阶跃信号时系统响应曲线。state.m程序以下：
p=I*(M+m)+M*m*L^2;
>> 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]

A =

0 1.0000 0 0
0 -0.0883 0.6293 0
0 0 0 1.0000
0 -0.2357 27.8285 0

>> B=[0;
(I+m*L^2)/p;
0;
m*L/p]

B =

0
0.8832
0
2.3566

>> C=[1 0 0 0;
0 0 1 0]
D=[0;
0]

C =

1 0 0 0
0 0 1 0

D =

0
0
matlab仿真开环阶跃响应曲线以下图所表示，系统并不稳定。

图2.2 系统开环阶跃响应曲线
三、 倒立摆系统PID控制算法设计

1. 试验要求和目标
l 要求：
设计PID控制器，使适当在小车上施加1N脉冲信号时，闭环系统响应指标为：
（1）稳定时间小于5秒
（2）稳态时摆杆和垂直方向夹角改变小于0.1 弧度
并作PID控制算法MATLAB仿真
l 目标：深入熟悉PID控制器设计方法，步骤，和P、I、D三参数调整方法。
2. 理论分析
l PID控制原理
在模拟控制系统中，控制器最常见控制规律是PID控制。常规PID控制系统原理框图以下图所表示。系统由模拟PID控制器KD(S)和被控对象G(S)组成。

PID控制器是一个线性控制器，它依据给定值和实际输出值组成控制偏差
将偏差百分比(P)、积分(I)和微分(D)经过线性组合组成控制量，对被控对象进行控制，故称PID控制器。其控制规律为

或写成传输函数形式

式中：——百分比系数；——积分时间常数；——微分时间常数。
在控制系统设计和仿真中，也将传输函数写成

式中：——百分比系数；——积分系数；——微分系数。
简单说来，PID控制器各校正步骤作用以下：
（1）百分比步骤：成百分比地反应控制系统偏差信号，偏差一旦产生，控制器立即产生控制作用，以降低偏差。
（2）积分步骤：关键用于消除稳态误差，提升系统型别。积分作用强弱取决于积分时间常数，越大，积分作用越弱，反之则越强。
（3）微分步骤：反应偏差信号改变趋势（改变速率），并能在偏差信号值变得太大之前，在系统中引入一个有效早期修正信号，从而加紧系统动作速度，减小调整时间。
l 摆杆角度控制
这个控制问题和我们以前碰到标准控制问题有些不一样，在这里输出量为摆杆位置，它初始位置为垂直向上，我们给系统施加一个扰动，观察摆杆响应。系统框图以下：

图中是控制器传输函数，是被控对象传输函数。
考虑到输入，结构图能够很轻易地变换成

该系统输出为

其中： ——被控对象传输函数分子项
——被控对象传输函数分母项
——PID控制器传输函数分子项
——PID控制器传输函数分母项
被控对象传输函数是

其中
PID控制器传输函数为

只需调整PID控制器参数，就能够得到满意控制效果。
l 小车位置控制
小车位置作为输出时，系统框图以下：

其中，是摆杆传输函数，是小车传输函数。
因为输入信号，所以能够把结构图转换成：

其中，反馈环代表我们前面设计摆杆控制器。
以后框图我们能够看出此处只对摆杆角度进行了控制，并没有对小车位置进行控制。
小车位置输出为：

其中，，，，分别代表被控对象1和被控对象2传输函数分子和分母。和代表PID控制器传输函数分子和分母。下面我们来求，依据前面试验二推导，有

能够推出小车位置传输函数为

其中
能够看出， ==，小车算式能够简化成：

3. PID控制算法MATLAB仿真
实际系统参数以下：
M 小车质量 1.096 Kg
m 摆杆质量 0.109 Kg
b 小车摩擦系数 0 .1N/m/sec
l 摆杆转动轴心到杆质心长度 0.25m
I 摆杆惯量 0.0034 kg*m*m
F 加在小车上力
x 小车位置
T 采样时间
摆杆matlab仿真程序代码以下：

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;
num1=[m*L/q 0 0];
den1=[1 b*(I+m*L^2)/q -(M+m)*m*g*L/q -b*m*g*L/q 0];
Kp=1;
Ki=1;
Kd=1;
numPID=[ Kd Kp Ki];
denPID=[1 0];
num=conv(num1,denPID);
den=polyadd(conv(denPID,den1),conv(numPID,num1));
[r,p,k]=residue(num,den);
s=p
t=0:0.005:5;
impulse(num,den,t)
axis([0 2 0 10])

运行程序得到：
s =

-6.4161
3.9693
0.0019
0
0
并得到仿真图像以下：

图3.1 kp=ki=kd=1时仿真响应图
可见此时系统并不稳定，此时应该首先调整kp，观察其响应改变：
讲kp设置为150，得到并观察响应图以下：

s =

-1.2224 +18.0044i
-1.2224 -18.0044i
-0.0000
-0.0000
-0.0000

图3-2，kp=150系统仿真图
可见此时系统两个闭环极点均在S平面做平面，系统稳定，系统稳定时间约为4秒，满足要求。此时系统有极小静态误差，依据系统对于精度要求可酌情考虑是否添加积分控制，本文添加积分控制。
将积分参数设为5，得到并观察闭环响应图。
在笔者经过数次尝试以后，发觉积分控制对于系统响应调整作用极小，笔者给出当积分参数分别设为10和50响应图以下：

图3-3，ki=10响应

图3-4，ki=50系统响应
积分作用通常是用来调整系统静态误差，使之达成需要范围，不过此处显著积分作用对系统影响不大，并了解到被控对象特征属于改变快类型，应该考虑改变微分控制，即使微分控制在实际系统中利用并不多见。
笔者将微分作用参数设置为10,20,50观察其效果图。

图3-5，kp=150,ki=50,kd=10仿真图像

图3-6，kp=150,ki=50,kd=20仿真图像

图3-7，kp=150,ki=50,kd=50仿真图像
当微分效果加上去时候，系统闭环仿真图像结果得到了质改善，瞬间替换了超调，不稳定，响应时间也快速降到了0.5秒，稳定时间在1秒，完美地完成了任务。其效果已经不能简单用好来形容，不过微分作用并非如此普及，而且每次全部效果如此良好，要依据不一样对象来判定用什么作用。必需要说是，微分作用在物理实现中是并不轻易，假如只有百分比调整和积分调整就能达成预想效果，那就不要使用微分调整。

4. 小车位置控制算法仿真

pid2.m是仿真小车位置改变m文件，文件以下：
% 小车位置PID控制
% 输入倒立摆传输函数 G1(s)=num1/den1,G2(s)=num2/den2
M = 1.096; m = 0.109; b = 0.1; I = 0.0034; g = 9.8; l= 0.25;
q = (M+m)*(I+m*l^2) -(m*l)^2;
num1 = [m*l/q 0 0];
den1 = [1 b*(I+m*l^2)/q -(M+m)*m*g*l/q -b*m*g*l/q 0];
num2 = [-(I+m*l^2)/q 0 m*g*l/q];
den2 = den1;
% 输入控制器PID数学模型 Gc(s)=numPID/denPID
Kp = 150; Ki = 50; Kd = 50;
numPID = [Kd Kp Ki];
denPID = [1 0];
% 计算闭环系统传输函数G(s)=num/den
% 多项式相乘
num = conv(num2,denPID);
% 多项式相加
den = polyadd(conv(denPID,den2),conv(numPID,num1 ));
% 求闭环系统极点
[r,p,k] = residue(num,den);
% 显示闭环系统极点
s = p
% 求取多项式传函脉冲响应
t=0:0.005:5;
impulse(num,den,t)
% 显示范围：横坐标0-5，纵坐标0-10，此条语句参数可依据仿真输出曲线调整
axis([0 5 -0.1 0.5])
grid

此时系统取Kp=150，Ki=50，Kd=50，阶跃响应仿真曲线以下图所表示：
s =

-115.0953
-2.4030
-0.4177
0
0

图3-8，小车位置仿真
由仿真结果能够看出，当摆杆角度处于很好闭环控制下时，小车位置即使处于失控状态，不过上升速度不快。

四、倒立摆系统最优控制算法设计
1. 设计目标和要求
现代控制理论最突出特点就是将控制对象用状态空间表示式形式表示出来，这么便于对多输入多输出系统进行分析和设计。线性二次型最优控制算法（）是现代控制理论中一个关键、基础方法，算法目标是在一定性能指标下，使系统控制效果最好，即利用最少控制能量，来达成最小状态误差。本章关键利用最优控制算法实现对一阶倒立摆系统摆杆角度和小车位置同时控制。
设计目标：
学习怎样使用状态空间法设计系统控制算法。
设计要求：
用状态空间法设计控制器，使适当在小车上施加0.2N阶跃信号时，闭环系统响应 指标为：
（1）摆杆角度和小车位移稳定时间小于5秒
（2）上升时间小于1秒
（3）超调量小于20度（0.35弧度）
（4）稳态误差小于2%。

2. 最优控制器设计
在PID调整中，我们输入是脉冲量，而且在设计控制器时，只对摆杆角度进行控制，而不考虑小车位移。然而，对一个倒立摆系统来说，把它作为单输出系统是不符合实际，假如把系统看成多输出系统话，用状态空间法分析要相对简单部分，在这一章我们将设计一个对摆杆位置和小车位移全部进行控制控制器。
系统状态方程为

在倒立摆相关参数为：
小车质量 1.096Kg
摆杆质量 0.2109Kg
小车摩擦系数 0.1 N/m/sec
摆杆转动轴心到杆质心长度 0.25 m
摆杆惯量 0.0034 kg*m*m
采样时间 0.005秒
条件下，状态方程系数矩阵以下：
；；；
最优控制前提条件是系统是能控，下面来判定一下系统能控能观性。
l Matlab仿真程序以下：
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]
Qc=ctrb(A,B);//判定能控性
K=rank(Qc)
Qo=obsv(A,C)；//判定能观性
I=rank(Qo)
1. Matlab仿真结果为：

K =

4
I =

4
即：系统能控矩阵秩 。
系统能观矩阵秩 。
故系统是能控能观。所以能够给系统加上最优控制器使得系统闭环稳定，且满足暂态性能指标。在利用线性二次型最优控制算法进行控制器设计时，关键目标就是取得反馈向量值。由上一小节推导知道，设计系统状态反馈控制器时，一个关键问题就是二次型性能指标泛函中加权矩阵和选择。为了使问题简化及加权矩阵含有比较明确物理意义，我们将取为对角阵。假设
；
这么得到性能指标泛函为

由上式能够看出，是正确平方加权，相对增加就意味着正确要求相对其它状态变量严格，在性能指标中比重大，偏差状态相对减小。是对控制量平方加权，当相对较大时，意味着控制费用增加，使得控制能量较小，反馈减弱，而取值较小时，系统控制费用减小，反馈增加，系统动态响应快速。
考虑到一阶倒立摆系统在运行过程中，关键被控量为系统输出量和，所以在选择加权对角阵各元素值时，因为代表小车位置权重，而是摆杆角度权重，所以只选择、，而。
选择和时需要注意多个方面：
（1）因为我们采取系统模型是线性化结果，为使系统个状态量能够在线性范围工作，要求各状态量不应过大。
（2）闭环系统最好能有一对共轭复数极点，这么有利于克服系统非线性摩擦，但系统主导极点模不应太大以免系统频带过宽，使得系统对噪声太敏感，以致系统不能正常工作。
（3）加权矩阵减小，会造成大控制能量，应注意控制大小，不要超出系统实施机构能力，使得放大器处于饱和状态。
控制系统以下图所表示，图中R是施加在小车上阶跃输入，四个状态量、、和分别代表小车位移、小车速度、摆杆位置和摆杆角速度，输出包含小车位置和摆杆角度。我们要设计一个控制器，使适当给系统施加一个阶跃输入时，摆杆会摆动，然后仍然回到垂直位置，小车抵达新命令位置。

2） 系统仿真
M = 0.5; m = 0.2; b = 0.1; I = 0.006; g = 9.8; l = 0.3;
p = I*(M+m)+M*m*l^2;
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]; p = eig(A);
% 求向量K
x = 5000; y = 100;
Q = [x 0 0 0;
0 0 0 0;
0 0 y 0
0 0 0 0];
R = 1; K = lqr(A,B,Q,R)
% 计算LQR控制矩阵
Ac = [(A-B*K)]; Bc = [B]; Cc = [C]; Dc = [D];
% 计算增益Nbar
Cn = [1 0 0 0];
Nbar = rscale(A,B,Cn,0,K);
Bcn = [Nbar*B];
% 求阶跃响应并显示，小车位置为虚线，摆杆角度为实线
T = 0:0.005:5;
U = 0.2*ones(size(T));
[Y,X] = Lsim(Ac,Bcn,Cc,Dc,U,T);
plot(T,Y(:,1),':',T,Y(:,2),'-')
legend('Cart Position','Pendulum Angle')
grid
文件中用到求取输入输出匹配系数函数rscale，它不是Matlab 工具，所以必需把它拷贝到rscale.m文件中, 并把该文件和源文件一起拷贝到MATLAB工作区。rscale.m文件以下：
% 求取输入输出匹配系数
function[Nbar] = rscale(A,B,C,D,K)
s = size(A,1);
Z = [zeros([1,s]) 1];
N = inv([A,B;C,D])*Z';
Nx = N(1:s);
Nu = N(1+s);
Nbar = Nu + K*Nx;
用函数来计算，运行程序，得到：
K =

-70.7107 -40.6531 125.7702 24.3770

图4-1，系统仿真图

即，能够看出，实际上和向量中和小车位置对应那一项相等。
此时系统响应曲线以下，小车位置跟踪输入信号；而且，摆杆超调足够小，稳态误差满足要求，上升时间和稳定时间也符合设计指标。

五、 参考文件
倒立摆课程设计指导书 曲延滨 3月 哈尔滨工业大学出版社