课程设计说明书
课程名称： 控制系统课程设计 设计题目：一阶倒立摆控制器设计 院 系： 信息与电气工程学院 班 级： 设 计 者： 学 号： 指导教师： 设计时间：2013年2月25日到2013年3月8号
课程设计（论文）任务书
指导教师签字： 系（教研室）主任签字：
2013年 3月 5日
目录
一、 建立一阶倒立摆数学模型 . ............................................................................................ 4
1. 一阶倒立摆的微分方程模型 . ..................................................................................... 4 2. 一阶倒立摆的传递函数模型 . ..................................................................................... 6 3. 一阶倒立摆的状态空间模型 . ..................................................................................... 7 二、 一阶倒立摆matlab 仿真 ................................................................................................ 9 三、 倒立摆系统的PID 控制算法设计 .............................................................................. 13 四、倒立摆系统的最优控制算法设计 . ................................................................................ 23 五、 总结 . .............................................................................................. 错误！未定义书签。 六、 参考文献 . ...................................................................................................................... 29
一、建立一阶倒立摆数学模型
首先建立一阶倒立摆的物理模型。在忽略空气阻力和各种摩擦之后， 可将直线一级倒立摆系统抽象成小车和匀质杆组成的系统，如图1所示。
系统内部各相关参数定义如下： M 小车质量 m 摆杆质量 b 小车摩擦系数
l 摆杆转动轴心到杆质心的长度 I 摆杆惯量
F 加在小车上的力 x 小车位置
φ 摆杆与垂直向上方向的夹角 θ 摆杆与垂直向下方向的夹角（考虑到摆杆初始位置为竖直向下）
1. 一阶倒立摆的微分方程模型
对一阶倒立摆系统中的小车和摆杆进行受力分析，其中，N 和 P为小车与摆杆相互作用力的水平和垂直方向的分量。
图 1-2 小车及摆杆受力图
分析小车水平方向所受的合力，可以得到以下方程：
（1-1）
由摆杆水平方向的受力进行分析可以得到下面等式：
（1-2）
即：
把这个等式代入式(1-1)中，就得到系统的第一个运动方程：
（1-4）
为了推出系统的第二个运动方程，我们对摆杆垂直方向上的合力进行分析，可以得到下面方程：
（1-3）
（1-5）
即：
（1-6）
力矩平衡方程如下：
（1-7）
由于面有负号。
所以等式前
合并这两个方程，约去 P和 N，得到第二个运动方程：
（1-8）
设
，（φ是摆杆与垂直向上方向之间的夹角），假设φ
d θ2
) =0。用u 代表被控对dt
象的输入力F ，利用上述近似进行线性化得直线一阶倒立摆的微分方程为：
弧度， 则可以进行近似处理：cos θ=-1, sin θ=-φ, (
（1-9）
2. 一阶倒立摆的传递函数模型
对式(1-9)进行拉普拉斯变换，得：
（2-1）
注意：推导传递函数时假设初始条件为 0。
由于输出为角度φ，求解方程组的第一个方程，可得：
（2-2）
或
（2-3）
如果令错误！未找到引用源。，则有：
（2-4）
把上式代入方程组（2-1）的第二个方程，得：
（2-5）
整理后得到传递函数：
（2-6）
其中
。
3. 一阶倒立摆的状态空间模型
设系统状态空间方程为：
（3-1）
方程组（2-9）对x 解代数方程，得到解如下：
..
（3-1）
整理后得到系统状态空间方程：
（3-2）
（3-3）
1
摆杆的惯量为I ml 2，代入（1-9）的第一个方程为：
3
得：
化简得：
设x =[x x φ
.
.
（3-4）
φ],u =x 则有：
'
..
（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-5）
二、一阶倒立摆matlab 仿真
实际系统参数如下，按照上面给出的例子求系统的传递函数、状态空间方程，并进行脉冲响应和阶跃响应的matlab 仿真。
M 小车质量 摆杆质量
1.096Kg 0.109Kg
0.1 N/m/sec
m
b
小车摩擦系数
l 摆杆转动轴心到杆质心的长度 0.25m
I 摆杆惯量 0.0034 kg*m*m
T 采样时间 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
由此可知，系统传递函数的多项式表达式为：
Φ(s ) 2.3566s 2
（2-1） G (s )==4
32
U (s ) s +0.0883s -27.8285s -2.3094s
系统的开环极点为（s):s 1=-5. 2780、s 2=5. 2727、s 3=-0. 0830、s 4=0，由于有一个开环极点位于S 平面的右半部，开环系统并不是稳定的。 系统的脉冲响应如下，由图也可见，系统并不稳定。
Impulse Response
A m p l i t u d e
0.5
Time (sec)
11.5
图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. 实验要求与目的 ● 要求：
设计PID 控制器，使得当在小车上施加1N 的脉冲信号时，闭环系统的响应指标为：
（1）稳定时间小于5秒
（2）稳态时摆杆与垂直方向的夹角变化小于0.1 弧度 并作PID 控制算法的MATLAB 仿真
● 目的：进一步熟悉PID 控制器的设计方法，步骤，以及P 、I 、D 三参数的
调节方法。 2. 理论分析
● PID 控制原理
在模拟控制系统中，控制器最常用的控制规律是PID 控制。常规PID 控制系统原理框图如下图所示。系统由模拟PID 控制器KD(S)和被控对象G(S)组成。
PID 控制器是一种线性控制器，它根据给定值r (t ) 与实际输出值y (t ) 构成控制偏差e (t )
e (t ) =r (t ) -y (t )
将偏差的比例(P)、积分(I)和微分(D)通过线性组合构成控制量，对被控对象进行控制，故称PID 控制器。其控制规律为
⎡1
u (t ) =K P ⎢e (t ) +
T I
⎣
t
⎰e (t ) dt +T D
de (t ) ⎤
⎥ dt ⎦
或写成传递函数的形式
⎛⎫U (s ) 1 G (s ) ==K P 1++T D s ⎪⎪ E (s ) T s I ⎝⎭
式中：K P ——比例系数；T I ——积分时间常数；T D ——微分时间常数。 在控制系统设计和仿真中，也将传递函数写成
G (s ) =
K U (s )
=K P +I +K D s E (s ) s
式中：K P ——比例系数；K I ——积分系数；K D ——微分系数。 简单说来，PID 控制器各校正环节的作用如下：
（1）比例环节：成比例地反映控制系统的偏差信号e (t ) ，偏差一旦产生，控制器立即产生控制作用，以减少偏差。
（2）积分环节：主要用于消除稳态误差，提高系统的型别。积分作用的强弱取决于积分时间常数T I ，T I 越大，积分作用越弱，反之则越强。
（3）微分环节：反映偏差信号的变化趋势（变化速率），并能在偏差信号值变得太大之前，在系统中引入一个有效的早期修正信号，从而加快系统的动作速度，减小调节时间。 摆杆角度控制
这个控制问题和我们以前遇到的标准控制问题有些不同，在这里输出量为摆杆的位置，它的初始位置为垂直向上，我们给系统施加一个扰动，观察摆杆的响应。系统框图如下：
图中KD (s ) 是控制器传递函数，G (s ) 是被控对象传递函数。 考虑到输入r (s ) =0，结构图可以很容易地变换成
该系统的输出为
num
G (s ) den y (s ) =F (s ) =F (s )
(numPID )(num ) 1+KD (s ) G (s )
1+
(denPID )(den )
num (denPID )
=F (s ) (denPID )(den ) +(numPID )(num )
其中： num ——被控对象传递函数的分子项
den ——被控对象传递函数的分母项
numPID ——PID 控制器传递函数的分子项 denPID ——PID 控制器传递函数的分母项
被控对象的传递函数是
ml 2
s
Φ(s ) q num
==U (s ) den b (I +ml 2) 3(M +m ) mgl 2bmgl 4
s +s -s -s
q q q
其中 q =[(M +m )(I +ml ) -(ml ) ] PID 控制器的传递函数为
2
2
K I K D s 2+K P s +K I numPID
KD (s ) =K D s +K P +==
s s denPID
只需调节PID 控制器的参数，就可以得到满意的控制效果。 小车位置控制
小车位置作为输出时，系统框图如下：
其中，G 1(s ) 是摆杆传递函数，G 2(s ) 是小车传递函数。
由于输入信号r (s ) =0，所以可以把结构图转换成：
其中，反馈环代表我们前面设计的摆杆的控制器。
从此框图我们可以看出此处只对摆杆角度进行了控制，并没有对小车位置进行控制。 小车位置输出为：
num 2
G 2(s ) den 2
X (s ) =F (s ) =F (s )
(numPID )(num ) 1+KD (s ) G 1(s ) 11+
(denPID )(den 1) =
(num 2)(denPID )(den 1)
F (s )
(denPID )(den 1)(den 2) +(numPID )(num 1)(den 2)
其中，num 1，den 1，num 2，den 2分别代表被控对象1和被控对象2传递函数的分子和分母。numPID 和denPID 代表PID 控制器传递函数的分子和分母。下面我们来求G 2(s ) ，根据前面实验二的推导，有
(I +ml 2) g
X (s ) =[-2]Φ(s )
ml s
可以推出小车位置的传递函数为
(I +ml 2) 2mgl
s -
X (s ) q q G 2(s ) ==2U (s ) b (I +ml ) 3(M +m ) mgl 2bmgl
s 4+s -s -s
q q q
其中 q =[(M +m )(I +ml ) -(ml ) ]
可以看出，den 1 =den 2=den ，小车的算式可以简化成：
2
2
X (s ) =
(num 2)(denPID )
F (s )
(denPID )(den ) +k (numPID )(num 1)
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
并得到仿真图像如下：
Impulse Response
109
876A m p l i t u d e
543210
00.20.40.60.81Time (sec)
1.21.41.61.82
图3.1 kp=ki=kd=1时的仿真响应图
可见此时系统并不稳定，此时应该首先调整kp ，观察其响应的变化： 讲kp 设置为150，得到并观察响应图如下： s =
-1.2224 +18.0044i -1.2224 -18.0044i -0.0000 -0.0000 -0.0000
Impulse Response
0.20.15
0.10.05
A m p l i t u d e
-0.05-0.1
-0.15
-0.2
00.511.52Time (sec)
2.533.54
图3-2，kp=150系统仿真图
可见此时系统两个闭环极点均在S 平面做平面，系统稳定，系统稳定时间约为4秒，满足要求。此时系统有极小的静态误差，根据系统对于精度的要求可酌情考虑是否添加积分控制，本文添加积分控制。 将积分参数设为5，得到并观察闭环响应图。
在笔者经过多次尝试之后，发现积分控制对于系统响应的调节作用极小，笔者给出当积分参数分别设为10和50的响应图如下：
Impulse Response
0.20.15
0.10.05A m p l i t u d e
-0.05-0.1
-0.15
-0.2
00.511.522.5Time (sec)
33.544.55
图3-3，ki=10的响应
Impulse Response
0.20.15
0.10.05A m p l i t u d e
-0.05-0.1
-0.15
-0.2
00.511.522.5Time (sec)
33.544.55
图3-4，ki=50系统的响应
积分作用通常是用来调整系统的静态误差，使之达到需要的范围，但是此处明显积分作用对系统的影响不大，并了解到被控对象的特性属于变化快的类型，应该考虑改变微分控制，虽然微分控制在实际系统中运用并不多见。 笔者将微分作用参数设置为10,20,50观察其效果图。
Impulse Response
0.20.15
0.10.05A m p l i t u d e
-0.05-0.1
-0.15
-0.2
00.511.522.5Time (sec)
33.544.55
图3-5，kp=150,ki=50,kd=10的仿真图像
Impulse Response
0.20.15
0.10.05
A m p l i t u d e
-0.05-0.1
-0.15
-0.2
00.511.522.5Time (sec)
33.544.55
图3-6，kp=150,ki=50,kd=20的仿真图像
Impulse Response
0.20.15
0.10.05A m p l i t u d e
-0.05-0.1
-0.15
-0.2
00.511.522.5Time (sec)
33.544.55
图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; 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 =
I = 0.0034; g = 9.8; l= 0.25;
-115.0953 -2.4030 -0.4177 0
Impulse Response
0.4
0.3
A m p l i t u d e
0.2
0.1
-0.1
00.511.522.5Time (sec)
33.544.55
图3-8，小车位置仿真
由仿真结果能够看出，当摆杆角度处于很好的闭环控制下时，小车位置虽然处于失控状态，但是上升速度不快。
四、倒立摆系统的最优控制算法设计
1. 设计目的与要求
现代控制理论的最突出特点就是将控制对象用状态空间表达式的形式表示出来，这样便于对多输入多输出系统进行分析和设计。线性二次型最优控制算法（LQR ）是现代控制理论中一种重要的、基本的方法，LQR 算法的目的是在一定的性能指标下，使系统的控制效果最佳，即利用最少的控制能量，来达到最小的状态误差。本章主要利用最优控制算法实现对一阶倒立摆系统的摆杆角度和小车位置的同时控制。 设计目的：
学习如何使用状态空间法设计系统的控制算法。 设计要求：
用状态空间法设计控制器，使得当在小车上施加0.2N 的阶跃信号时，闭环系统的响应 指标为：
（1）摆杆角度θ和小车位移x 的稳定时间小于5秒 （2）x 的上升时间小于1秒
（3）θ的超调量小于20度（0.35弧度） （4）稳态误差小于2%。
2. 最优控制器的设计
在PID 调节中，我们的输入是脉冲量，并且在设计控制器时，只对摆杆角度进行控制，而不考虑小车的位移。然而，对一个倒立摆系统来说，把它作为单输出系统是不符合实际的，如果把系统当作多输出系统的话，用状态空间法分析要相对简单一些，在这一章我们将设计一个对摆杆位置和小车位移都进行控制的控制器。 系统状态方程为
=AX +Bu X
Y =CX +Du
在倒立摆相关参数为：
M 小车质量 1.096Kg
摆杆质量 m
b 小车摩擦系数
0.2109Kg
0.1 N/m/sec
l 摆杆转动轴心到杆质心的长度 0.25 m
I 摆杆惯量 T 采样时间
0.0034 kg*m*m 0.005秒
的条件下，状态方程系数矩阵如下：
1⎡0
⎢00. 0883A =⎢
⎢00⎢
⎣0-0. 2357
⎤⎡0⎤
⎢0. 8832⎥⎡1000⎤01. 0000⎥⎥；B =⎢⎥；C =⎢⎥；0010⎢0⎥01. 0000⎥⎣⎦
⎥⎢⎥
27. 82850⎦⎣2. 3566⎦0
⎡0⎤
D =⎢⎥
⎣0⎦
最优控制的前提条件是系统是能控的，下面来判断一下系统的能控能观性。 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
即：系统的能控矩阵的秩 rank [B AB A 2B A 3B ]=4。 系统的能观矩阵的秩 rank [C CA CA 2CA 3]=4。
故系统是能控能观的。因此可以给系统加上最优控制器使得系统闭环稳定，且满足暂态性能指标。在运用线性二次型最优控制算法进行控制器设计时，主要的目的就是获得反馈向量K 的值。由上一小节的推导知道，设计系统状态反馈控制器时，一个关键的问题就是二次型性能指标泛函中加权矩阵Q 和R 的选取。为了使问题简化及加权矩阵具有比较明确的物理意义，我们将Q 取为对角阵。假设
⎡Q 110
⎢0Q
22
Q =⎢
⎢00⎢
0⎣0
00Q 330
0⎤
0⎥⎥；R =[r ] 0⎥⎥Q 44⎦
这样得到的性能指标泛函为
∞
J =⎰Q 11x 1+Q 22x 2+Q 33x 3+Q 44x 4+ru 2dt
[
2222
]
由上式可以看出，Q ii 是对x i 的平方的加权，Q ii 的相对增加就意味着对x i 的要求相对其它状态变量严格，在性能指标中的比重大，x i 的偏差状态相对减小。r 是对控制量u 的平方加权，当r 相对较大时，意味着控制费用增加，使得控制能量较小，反馈减弱，而r 取值较小时，系统控制费用减小，反馈增加，系统动态响应迅速。
考虑到一阶倒立摆系统在运行过程中，主要的被控量为系统的输出量x 和φ，因此在选取加权对角阵Q 的各元素值时，由于Q 11代表小车位置的权重，而Q 33是摆杆角度的权重，所以只选取Q 11、Q 33，而Q 22=Q 44=0。 选取Q 和R 时需要注意的几个方面：
（1）由于我们采用的系统模型是线性化的结果，为使系统个状态量能够在线性范围工作，要求各状态量不应过大。
（2）闭环系统最好能有一对共轭复数极点，这样有利于克服系统的非线性摩擦，但系统主导极点的模不应太大以免系统频带过宽，使得系统对噪声太敏感，以致系统不能正常工作。
（3）加权矩阵R 的减小，会导致大的控制能量，应注意控制U 的大小，不要超过系统执行机构的能力，使得放大器处于饱和状态。
、控制系统如下图所示，图中R 是施加在小车上的阶跃输入，四个状态量x 、x 'φ和φ 分别代表小车位移、小车速度、摆杆位置和摆杆角速度，输出y =[x , φ]包
括小车位置和摆杆角度。我们要设计一个控制器，使得当给系统施加一个阶跃输入时，摆杆会摆动，然后仍然回到垂直位置，小车到达新的命令位置。
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;
用函数rscale 来计算Nbar ，运行程序，得到： K =
-70.7107 -40.6531 125.7702 24.3770
0.25
0.20.150.10.050-0.05-0.1-0.150
0.511.522.533.544.55
图4-1，系统仿真图
即Nbar =rscale (A , B , Cn , 0, K ) =-70. 7107，可以看出，实际上Nbar 和K 向量中与小车位置x 对应的那一项相等。
此时系统的响应曲线如下，小车位置跟踪输入信号；并且，摆杆超调足够小，稳态误差满足要求，上升时间和稳定时间也符合设计指标。
五、参考文献
倒立摆课程设计指导书 曲延滨 2011年3月 哈尔滨工业大学出版社
课程设计说明书
课程名称： 控制系统课程设计 设计题目：一阶倒立摆控制器设计 院 系： 信息与电气工程学院 班 级： 设 计 者： 学 号： 指导教师： 设计时间：2013年2月25日到2013年3月8号
课程设计（论文）任务书
指导教师签字： 系（教研室）主任签字：
2013年 3月 5日
目录
一、 建立一阶倒立摆数学模型 . ............................................................................................ 4
1. 一阶倒立摆的微分方程模型 . ..................................................................................... 4 2. 一阶倒立摆的传递函数模型 . ..................................................................................... 6 3. 一阶倒立摆的状态空间模型 . ..................................................................................... 7 二、 一阶倒立摆matlab 仿真 ................................................................................................ 9 三、 倒立摆系统的PID 控制算法设计 .............................................................................. 13 四、倒立摆系统的最优控制算法设计 . ................................................................................ 23 五、 总结 . .............................................................................................. 错误！未定义书签。 六、 参考文献 . ...................................................................................................................... 29
一、建立一阶倒立摆数学模型
首先建立一阶倒立摆的物理模型。在忽略空气阻力和各种摩擦之后， 可将直线一级倒立摆系统抽象成小车和匀质杆组成的系统，如图1所示。
系统内部各相关参数定义如下： M 小车质量 m 摆杆质量 b 小车摩擦系数
l 摆杆转动轴心到杆质心的长度 I 摆杆惯量
F 加在小车上的力 x 小车位置
φ 摆杆与垂直向上方向的夹角 θ 摆杆与垂直向下方向的夹角（考虑到摆杆初始位置为竖直向下）
1. 一阶倒立摆的微分方程模型
对一阶倒立摆系统中的小车和摆杆进行受力分析，其中，N 和 P为小车与摆杆相互作用力的水平和垂直方向的分量。
图 1-2 小车及摆杆受力图
分析小车水平方向所受的合力，可以得到以下方程：
（1-1）
由摆杆水平方向的受力进行分析可以得到下面等式：
（1-2）
即：
把这个等式代入式(1-1)中，就得到系统的第一个运动方程：
（1-4）
为了推出系统的第二个运动方程，我们对摆杆垂直方向上的合力进行分析，可以得到下面方程：
（1-3）
（1-5）
即：
（1-6）
力矩平衡方程如下：
（1-7）
由于面有负号。
所以等式前
合并这两个方程，约去 P和 N，得到第二个运动方程：
（1-8）
设
，（φ是摆杆与垂直向上方向之间的夹角），假设φ
d θ2
) =0。用u 代表被控对dt
象的输入力F ，利用上述近似进行线性化得直线一阶倒立摆的微分方程为：
弧度， 则可以进行近似处理：cos θ=-1, sin θ=-φ, (
（1-9）
2. 一阶倒立摆的传递函数模型
对式(1-9)进行拉普拉斯变换，得：
（2-1）
注意：推导传递函数时假设初始条件为 0。
由于输出为角度φ，求解方程组的第一个方程，可得：
（2-2）
或
（2-3）
如果令错误！未找到引用源。，则有：
（2-4）
把上式代入方程组（2-1）的第二个方程，得：
（2-5）
整理后得到传递函数：
（2-6）
其中
。
3. 一阶倒立摆的状态空间模型
设系统状态空间方程为：
（3-1）
方程组（2-9）对x 解代数方程，得到解如下：
..
（3-1）
整理后得到系统状态空间方程：
（3-2）
（3-3）
1
摆杆的惯量为I ml 2，代入（1-9）的第一个方程为：
3
得：
化简得：
设x =[x x φ
.
.
（3-4）
φ],u =x 则有：
'
..
（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-5）
二、一阶倒立摆matlab 仿真
实际系统参数如下，按照上面给出的例子求系统的传递函数、状态空间方程，并进行脉冲响应和阶跃响应的matlab 仿真。
M 小车质量 摆杆质量
1.096Kg 0.109Kg
0.1 N/m/sec
m
b
小车摩擦系数
l 摆杆转动轴心到杆质心的长度 0.25m
I 摆杆惯量 0.0034 kg*m*m
T 采样时间 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
由此可知，系统传递函数的多项式表达式为：
Φ(s ) 2.3566s 2
（2-1） G (s )==4
32
U (s ) s +0.0883s -27.8285s -2.3094s
系统的开环极点为（s):s 1=-5. 2780、s 2=5. 2727、s 3=-0. 0830、s 4=0，由于有一个开环极点位于S 平面的右半部，开环系统并不是稳定的。 系统的脉冲响应如下，由图也可见，系统并不稳定。
Impulse Response
A m p l i t u d e
0.5
Time (sec)
11.5
图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. 实验要求与目的 ● 要求：
设计PID 控制器，使得当在小车上施加1N 的脉冲信号时，闭环系统的响应指标为：
（1）稳定时间小于5秒
（2）稳态时摆杆与垂直方向的夹角变化小于0.1 弧度 并作PID 控制算法的MATLAB 仿真
● 目的：进一步熟悉PID 控制器的设计方法，步骤，以及P 、I 、D 三参数的
调节方法。 2. 理论分析
● PID 控制原理
在模拟控制系统中，控制器最常用的控制规律是PID 控制。常规PID 控制系统原理框图如下图所示。系统由模拟PID 控制器KD(S)和被控对象G(S)组成。
PID 控制器是一种线性控制器，它根据给定值r (t ) 与实际输出值y (t ) 构成控制偏差e (t )
e (t ) =r (t ) -y (t )
将偏差的比例(P)、积分(I)和微分(D)通过线性组合构成控制量，对被控对象进行控制，故称PID 控制器。其控制规律为
⎡1
u (t ) =K P ⎢e (t ) +
T I
⎣
t
⎰e (t ) dt +T D
de (t ) ⎤
⎥ dt ⎦
或写成传递函数的形式
⎛⎫U (s ) 1 G (s ) ==K P 1++T D s ⎪⎪ E (s ) T s I ⎝⎭
式中：K P ——比例系数；T I ——积分时间常数；T D ——微分时间常数。 在控制系统设计和仿真中，也将传递函数写成
G (s ) =
K U (s )
=K P +I +K D s E (s ) s
式中：K P ——比例系数；K I ——积分系数；K D ——微分系数。 简单说来，PID 控制器各校正环节的作用如下：
（1）比例环节：成比例地反映控制系统的偏差信号e (t ) ，偏差一旦产生，控制器立即产生控制作用，以减少偏差。
（2）积分环节：主要用于消除稳态误差，提高系统的型别。积分作用的强弱取决于积分时间常数T I ，T I 越大，积分作用越弱，反之则越强。
（3）微分环节：反映偏差信号的变化趋势（变化速率），并能在偏差信号值变得太大之前，在系统中引入一个有效的早期修正信号，从而加快系统的动作速度，减小调节时间。 摆杆角度控制
这个控制问题和我们以前遇到的标准控制问题有些不同，在这里输出量为摆杆的位置，它的初始位置为垂直向上，我们给系统施加一个扰动，观察摆杆的响应。系统框图如下：
图中KD (s ) 是控制器传递函数，G (s ) 是被控对象传递函数。 考虑到输入r (s ) =0，结构图可以很容易地变换成
该系统的输出为
num
G (s ) den y (s ) =F (s ) =F (s )
(numPID )(num ) 1+KD (s ) G (s )
1+
(denPID )(den )
num (denPID )
=F (s ) (denPID )(den ) +(numPID )(num )
其中： num ——被控对象传递函数的分子项
den ——被控对象传递函数的分母项
numPID ——PID 控制器传递函数的分子项 denPID ——PID 控制器传递函数的分母项
被控对象的传递函数是
ml 2
s
Φ(s ) q num
==U (s ) den b (I +ml 2) 3(M +m ) mgl 2bmgl 4
s +s -s -s
q q q
其中 q =[(M +m )(I +ml ) -(ml ) ] PID 控制器的传递函数为
2
2
K I K D s 2+K P s +K I numPID
KD (s ) =K D s +K P +==
s s denPID
只需调节PID 控制器的参数，就可以得到满意的控制效果。 小车位置控制
小车位置作为输出时，系统框图如下：
其中，G 1(s ) 是摆杆传递函数，G 2(s ) 是小车传递函数。
由于输入信号r (s ) =0，所以可以把结构图转换成：
其中，反馈环代表我们前面设计的摆杆的控制器。
从此框图我们可以看出此处只对摆杆角度进行了控制，并没有对小车位置进行控制。 小车位置输出为：
num 2
G 2(s ) den 2
X (s ) =F (s ) =F (s )
(numPID )(num ) 1+KD (s ) G 1(s ) 11+
(denPID )(den 1) =
(num 2)(denPID )(den 1)
F (s )
(denPID )(den 1)(den 2) +(numPID )(num 1)(den 2)
其中，num 1，den 1，num 2，den 2分别代表被控对象1和被控对象2传递函数的分子和分母。numPID 和denPID 代表PID 控制器传递函数的分子和分母。下面我们来求G 2(s ) ，根据前面实验二的推导，有
(I +ml 2) g
X (s ) =[-2]Φ(s )
ml s
可以推出小车位置的传递函数为
(I +ml 2) 2mgl
s -
X (s ) q q G 2(s ) ==2U (s ) b (I +ml ) 3(M +m ) mgl 2bmgl
s 4+s -s -s
q q q
其中 q =[(M +m )(I +ml ) -(ml ) ]
可以看出，den 1 =den 2=den ，小车的算式可以简化成：
2
2
X (s ) =
(num 2)(denPID )
F (s )
(denPID )(den ) +k (numPID )(num 1)
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
并得到仿真图像如下：
Impulse Response
109
876A m p l i t u d e
543210
00.20.40.60.81Time (sec)
1.21.41.61.82
图3.1 kp=ki=kd=1时的仿真响应图
可见此时系统并不稳定，此时应该首先调整kp ，观察其响应的变化： 讲kp 设置为150，得到并观察响应图如下： s =
-1.2224 +18.0044i -1.2224 -18.0044i -0.0000 -0.0000 -0.0000
Impulse Response
0.20.15
0.10.05
A m p l i t u d e
-0.05-0.1
-0.15
-0.2
00.511.52Time (sec)
2.533.54
图3-2，kp=150系统仿真图
可见此时系统两个闭环极点均在S 平面做平面，系统稳定，系统稳定时间约为4秒，满足要求。此时系统有极小的静态误差，根据系统对于精度的要求可酌情考虑是否添加积分控制，本文添加积分控制。 将积分参数设为5，得到并观察闭环响应图。
在笔者经过多次尝试之后，发现积分控制对于系统响应的调节作用极小，笔者给出当积分参数分别设为10和50的响应图如下：
Impulse Response
0.20.15
0.10.05A m p l i t u d e
-0.05-0.1
-0.15
-0.2
00.511.522.5Time (sec)
33.544.55
图3-3，ki=10的响应
Impulse Response
0.20.15
0.10.05A m p l i t u d e
-0.05-0.1
-0.15
-0.2
00.511.522.5Time (sec)
33.544.55
图3-4，ki=50系统的响应
积分作用通常是用来调整系统的静态误差，使之达到需要的范围，但是此处明显积分作用对系统的影响不大，并了解到被控对象的特性属于变化快的类型，应该考虑改变微分控制，虽然微分控制在实际系统中运用并不多见。 笔者将微分作用参数设置为10,20,50观察其效果图。
Impulse Response
0.20.15
0.10.05A m p l i t u d e
-0.05-0.1
-0.15
-0.2
00.511.522.5Time (sec)
33.544.55
图3-5，kp=150,ki=50,kd=10的仿真图像
Impulse Response
0.20.15
0.10.05
A m p l i t u d e
-0.05-0.1
-0.15
-0.2
00.511.522.5Time (sec)
33.544.55
图3-6，kp=150,ki=50,kd=20的仿真图像
Impulse Response
0.20.15
0.10.05A m p l i t u d e
-0.05-0.1
-0.15
-0.2
00.511.522.5Time (sec)
33.544.55
图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; 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 =
I = 0.0034; g = 9.8; l= 0.25;
-115.0953 -2.4030 -0.4177 0
Impulse Response
0.4
0.3
A m p l i t u d e
0.2
0.1
-0.1
00.511.522.5Time (sec)
33.544.55
图3-8，小车位置仿真
由仿真结果能够看出，当摆杆角度处于很好的闭环控制下时，小车位置虽然处于失控状态，但是上升速度不快。
四、倒立摆系统的最优控制算法设计
1. 设计目的与要求
现代控制理论的最突出特点就是将控制对象用状态空间表达式的形式表示出来，这样便于对多输入多输出系统进行分析和设计。线性二次型最优控制算法（LQR ）是现代控制理论中一种重要的、基本的方法，LQR 算法的目的是在一定的性能指标下，使系统的控制效果最佳，即利用最少的控制能量，来达到最小的状态误差。本章主要利用最优控制算法实现对一阶倒立摆系统的摆杆角度和小车位置的同时控制。 设计目的：
学习如何使用状态空间法设计系统的控制算法。 设计要求：
用状态空间法设计控制器，使得当在小车上施加0.2N 的阶跃信号时，闭环系统的响应 指标为：
（1）摆杆角度θ和小车位移x 的稳定时间小于5秒 （2）x 的上升时间小于1秒
（3）θ的超调量小于20度（0.35弧度） （4）稳态误差小于2%。
2. 最优控制器的设计
在PID 调节中，我们的输入是脉冲量，并且在设计控制器时，只对摆杆角度进行控制，而不考虑小车的位移。然而，对一个倒立摆系统来说，把它作为单输出系统是不符合实际的，如果把系统当作多输出系统的话，用状态空间法分析要相对简单一些，在这一章我们将设计一个对摆杆位置和小车位移都进行控制的控制器。 系统状态方程为
=AX +Bu X
Y =CX +Du
在倒立摆相关参数为：
M 小车质量 1.096Kg
摆杆质量 m
b 小车摩擦系数
0.2109Kg
0.1 N/m/sec
l 摆杆转动轴心到杆质心的长度 0.25 m
I 摆杆惯量 T 采样时间
0.0034 kg*m*m 0.005秒
的条件下，状态方程系数矩阵如下：
1⎡0
⎢00. 0883A =⎢
⎢00⎢
⎣0-0. 2357
⎤⎡0⎤
⎢0. 8832⎥⎡1000⎤01. 0000⎥⎥；B =⎢⎥；C =⎢⎥；0010⎢0⎥01. 0000⎥⎣⎦
⎥⎢⎥
27. 82850⎦⎣2. 3566⎦0
⎡0⎤
D =⎢⎥
⎣0⎦
最优控制的前提条件是系统是能控的，下面来判断一下系统的能控能观性。 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
即：系统的能控矩阵的秩 rank [B AB A 2B A 3B ]=4。 系统的能观矩阵的秩 rank [C CA CA 2CA 3]=4。
故系统是能控能观的。因此可以给系统加上最优控制器使得系统闭环稳定，且满足暂态性能指标。在运用线性二次型最优控制算法进行控制器设计时，主要的目的就是获得反馈向量K 的值。由上一小节的推导知道，设计系统状态反馈控制器时，一个关键的问题就是二次型性能指标泛函中加权矩阵Q 和R 的选取。为了使问题简化及加权矩阵具有比较明确的物理意义，我们将Q 取为对角阵。假设
⎡Q 110
⎢0Q
22
Q =⎢
⎢00⎢
0⎣0
00Q 330
0⎤
0⎥⎥；R =[r ] 0⎥⎥Q 44⎦
这样得到的性能指标泛函为
∞
J =⎰Q 11x 1+Q 22x 2+Q 33x 3+Q 44x 4+ru 2dt
[
2222
]
由上式可以看出，Q ii 是对x i 的平方的加权，Q ii 的相对增加就意味着对x i 的要求相对其它状态变量严格，在性能指标中的比重大，x i 的偏差状态相对减小。r 是对控制量u 的平方加权，当r 相对较大时，意味着控制费用增加，使得控制能量较小，反馈减弱，而r 取值较小时，系统控制费用减小，反馈增加，系统动态响应迅速。
考虑到一阶倒立摆系统在运行过程中，主要的被控量为系统的输出量x 和φ，因此在选取加权对角阵Q 的各元素值时，由于Q 11代表小车位置的权重，而Q 33是摆杆角度的权重，所以只选取Q 11、Q 33，而Q 22=Q 44=0。 选取Q 和R 时需要注意的几个方面：
（1）由于我们采用的系统模型是线性化的结果，为使系统个状态量能够在线性范围工作，要求各状态量不应过大。
（2）闭环系统最好能有一对共轭复数极点，这样有利于克服系统的非线性摩擦，但系统主导极点的模不应太大以免系统频带过宽，使得系统对噪声太敏感，以致系统不能正常工作。
（3）加权矩阵R 的减小，会导致大的控制能量，应注意控制U 的大小，不要超过系统执行机构的能力，使得放大器处于饱和状态。
、控制系统如下图所示，图中R 是施加在小车上的阶跃输入，四个状态量x 、x 'φ和φ 分别代表小车位移、小车速度、摆杆位置和摆杆角速度，输出y =[x , φ]包
括小车位置和摆杆角度。我们要设计一个控制器，使得当给系统施加一个阶跃输入时，摆杆会摆动，然后仍然回到垂直位置，小车到达新的命令位置。
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;
用函数rscale 来计算Nbar ，运行程序，得到： K =
-70.7107 -40.6531 125.7702 24.3770
0.25
0.20.150.10.050-0.05-0.1-0.150
0.511.522.533.544.55
图4-1，系统仿真图
即Nbar =rscale (A , B , Cn , 0, K ) =-70. 7107，可以看出，实际上Nbar 和K 向量中与小车位置x 对应的那一项相等。
此时系统的响应曲线如下，小车位置跟踪输入信号；并且，摆杆超调足够小，稳态误差满足要求，上升时间和稳定时间也符合设计指标。
五、参考文献
倒立摆课程设计指导书 曲延滨 2011年3月 哈尔滨工业大学出版社