倒立摆实验报告(现代控制理论)

现代控制理论实验报告

——倒立摆

小组成员：

指导老师：

2013.5

实验一建立一级倒立摆的数学模型一、实验目的

学习建立一级倒立摆系统的数学模型，并进行Matlab仿真。二、实验内容

写出系统传递函数和状态空间方程，用Matlab进行仿真。三、Matlab源程序及程序运行的结果

（1）Matlab源程序见附页

（2）给出系统的传递函数和状态方程

（a）传递函数gs为摆杆的角度：

>> gs

Transfer function:

2.054 s

-----------------------------------

s^3 + 0.07391 s^2 - 29.23 s - 2.013

（b）传递函数gspo为小车的位移传递函数：

>> gspo

Transfer function:

0.7391 s^2 - 20.13

---------------------------------------

s^4 + 0.07391 s^3 - 29.23 s^2 - 2.013 s

（c）状态矩阵A,B,C,D：

>> sys

a =

x1 x2 x3 x4

x1 0 1 0 0

x2 0 -0.07391 0.7175 0

x3 0 0 0 1

x4 0 -0.2054 29.23 0

b =

u1

x1 0

x2 0.7391

x3 0

x4 2.054

c =

x1 x2 x3 x4

y1 1 0 0 0

y2 0 0 1 0

d =

u1

y1 0

y2 0

Continuous-time model.

（3）给出传递函数极点和系统状态矩阵A的特征值（a）传递函数gs的极点

>> P

P =

5.4042

-5.4093

-0.0689

（b）传递函数gspo的极点

>> Po

Po =

5.4042

-5.4093

-0.0689

（c）状态矩阵A的特征值

>> E

E =

-0.0689

5.4042

-5.4093

（4）给出系统开环脉冲响应和阶跃响应的曲线（a）开环脉冲响应曲线

（b）阶跃响应曲线

四、思考题

（1）由状态空间方程转化为传递函数，是否与直接计算传递函数相等？答：由状态空间方程转化为传递函数：

>> gso=tf(sys)

Transfer function from input to output...

0.7391 s^2 - 6.565e-016 s - 20.13

#1: ---------------------------------------

s^4 + 0.07391 s^3 - 29.23 s^2 - 2.013 s

2.054 s + 4.587e-016

#2: -----------------------------------

s^3 + 0.07391 s^2 - 29.23 s - 2.013

#1为gspo传递函数，#2为gs的传递函数

而直接得到的传递函数为：

>> gspo

Transfer function:

0.7391 s^2 - 20.13

---------------------------------------

s^4 + 0.07391 s^3 - 29.23 s^2 - 2.013 s

>> gs

Transfer function:

2.054 s

-----------------------------------

s^3 + 0.07391 s^2 - 29.23 s - 2.013

通过比较可以看到，gspo由状态空间方程转化的传递函数比直接得到的传递函数多了s的一次项，而6.565e-016非常小几乎可以忽略不计，因此可以认为两种方法得到的传递函数式相同的，同理传递函数gs也可以认为是相同的。（2）通过仿真表明开环系统是否稳定？请通过极点(特征值)理论来分析。答：开环系统不稳定

极点为：

>> P

P =

5.4042

-5.4093

-0.0689

>> Po

Po =

5.4042

-5.4093

-0.0689

由系统稳定性结论可知，极点若都分布在s平面的左半平面则系统稳定，而开环系统的极点有5.4042在右半平面。因此，开环系统不稳定。

（3）传递函数的极点和状态方程的特征值的个数、大小是否相等？如果不相等，请解释其原因。

传递函数gspo的极点和状态方程的特征值的个数、大小相等。但是传递函数gs的极点和状态方程的特征值个数不相等。因为存在零极点对消。

附录：（matlab程序）

clear all;

f1=0.001;

%实际系统参数

M=1.32;

m=0.132;

b=0.1;

l=0.27;

I=0.0032;

g=9.8;

T=0.02;

%求传递函数gs(输出为摆杆角度)和gspo(输出为小车位置)

q=(M+m)*(I+m*l^2)-(m*l)^2;

num=[m*l/q 0];

den=[1 b*(I+m*l^2)/q -(M+m)*m*g*l/q -b*m*g*l/q];

gs=tf(num,den);

numpo=[(I+m*l^2)/q 0 -m*g*l/q];

denpo=[1 b*(I+m*l^2)/q -(M+m)*m*g*l/q -b*m*g*l/q 0];

gspo=tf(numpo,denpo);

%求状态空间sys(A,B,C,D)

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*b*l/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];

sys=ss(A,B,C,D);

%通过传递函数求系统(摆杆角度和小车位置)的开环脉冲响应

t=0:T:5;

y1=impulse(gs,t);

y2=impulse(gspo,t);

figure(1);

plot(t,y2,'b',t,y1,'r');

xlabel('t/s');

ylabel('Position/m or Angle/rad');

axis([0 2 0 80]);

legend('Car Position','Pendulum Angle');

%将状态空间方程sys转化为传递函数gs0

gs0=tf(sys);

%通过状态方程求系统(摆杆角度和小车位置)的开环脉冲响应

t=0:T:5;

y=impulse(sys,t);

figure(2);

plot(t,y(:,1),t,y(:,2),'r');

xlabel('t/s');

ylabel('Position/m or Angle/rad');

axis([0 2 0 80]);

legend('Car Position','Pendulum Angle');

%通过传递函数求系统(摆杆角度和小车位置)的开环阶越响应

t=0:T:5;

y1=step(gs,t);

y2=step(gspo,t);

figure(3);

plot(t,y2,'b',t,y1,'r');

axis([0 2.5 0 80]);

xlabel('t/s');

ylabel('Position/m or Angle/rad');

legend('Car Position','Pendulum Angle');

%通过状态方程求系统(摆杆角度和小车位置)的开环阶越响应t=0:T:5;

y=step(sys,t);

figure(4);

plot(t,y(:,1),t,y(:,2),'r');

xlabel('t/s');

ylabel('Position/m or Angle/rad');

axis([0 2.5 0 80]);

legend('Car Position','Pendulum Angle');

%求传递函数极点

P=pole(gs);

Po=pole(gspo);

%求A的特征值

E=eig(A);

实验二倒立摆系统控制算法的状态空间法设计

一、实验目的

学习如何使用状态空间法设计系统的控制算法。

二、实验内容

用状态空间法设计控制器，使得当在小车上施加 0.2m 的阶跃信号时，闭环系统的响应指标为：

(1) 杆角度μ 和小车位移x的稳定时间小于5秒

(2) x的上升时间小于2秒2

(3) μ 的超调量小于20度（0.35弧度）

(4) 稳态误差小于4%.

三、Matlab源程序及程序执行结果

（1）Matlab源程序（见附录）

（2）程序执行结果

（a）k的值

>> K

K =

-14.1421 -12.1570 63.5837 11.8416

（b）反馈后的响应曲线

（3）给出无扰动时两次不同K 值下，小车的稳定位置P1和摆杆的稳定角度Pend1；

（a）>> K

K =

-14.142 -12.157 63.584 11.842

小车的稳定位置P1=-0.02

绿色的曲线为摆杆的稳定角度Pend1=0.001度(b)

>> K

K =

-14.1421 -12.1570 63.5837 11.8416

小车的稳定位置P1=-0.007

绿色的曲线为摆杆的稳定角度Pend1=0.0015度

（4）给出两次不同 K 值下，实际系统的响应曲线，并计算实验要求中的四项响应指标，并注意要利用实验三中统计出的响应时间延迟修正响应曲线。 ①K =

-14.1421 -12.1467 63.5825 11.8413

=(0.11-0.0825)/0.0925=29.7%

tp= (4100-3880)/1000*8.8=1.936s

tr=(4030-3880)/1000*8.8=1.32s ts=(4800-3880)/1000*8.8=8.096s

②K =

-14.1421 -12.1570 63.5837 11.8416

=(0.11-0.085)/0.092=27.17%

tp= (3025-2840)/1000*8.8=1.628s

tr=(2955-2840)/1000*8.8=1.012s ts=(4800-3900)/1000*8.8=7.92s

四、思考题

(1) 计算Ac 的特征值。 ①K =

-14.1421 -12.1467 63.5825 11.8413

②K =

-14.1421 -12.1570 63.5837 11.8416

(2) 通过仿真分析Q11和Q33的大小对控制效果的影响(Q11为Q阵的第(1; 1)个元素)：

?固定Q33 ，改变Q11

Q33=100 Q11=100（红）、500（蓝）、1000（绿）

从图中可以看出Q11增大,角度超调随着增大,位置的超调基本不变,但是响应时间缩短了。

?固定Q11 ，改变Q33

Q11=100 Q33=100（红）、1000（蓝）、2000（绿）

从图中可以看出Q33增大,角度超调减小,位置的超调基本不变,但是响应时间延长了。

附录：（matlab程序）

clear all;

f1=0.001;

%实际系统参数

% M=1.096;

% m=0.109;

% b=0.25;

% l=0.25;

% I=0.0034;

% g=9.8;

% T=0.001;

%求系统状态空间参数

M=1.32;

m=0.132;

b=0.22;

l=0.27;

I=0.0032;

g=9.8;

T=0.02;

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*b*l/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;

%求反馈向量K

R=1;

Q1=200;

Q2=0;

Q3=100;

Q=[Q1 0 0 0;0 Q2 0 0;0 0 Q3 0;0 0 0 0];

K=lqr(A,B,Q,R);

%求状态反馈后的系统sysstate

Ac=A-B*K;

Bc=B*K(1); %输入变换使输入与反馈的量纲匹配

sysstate=ss(Ac,Bc,C,D);

%对lqr控制系统进行仿真

t=0:T:5;

U=0.2*ones(size(t));

y=lsim(sysstate,U,t);

figure(1);

hold on;

plot(t,y(:,1),t,y(:,2),'r');

box on;

xlabel('t/s');

ylabel('Position/m or Angle/rad'); legend('Car Position','Pendulum Angle');

实验三研究倒立摆系统对信号的跟踪

一、实验目的

观察倒立摆对于不同输入信号的跟踪情况，加深对状态空间和状态反馈的理解。

二、实验内容

在平衡位置，分别设定下列三种信号，记录倒立摆的运动情况：

(1) 方波信号：频率0.2Hz，幅值0.05m

(2) 正弦波信号：频率0.2Hz，幅值0.05m

(3) 锯齿波信号：频率0.2Hz，幅值0.05m

三、Matlab源程序及程序执行结果

(1)Matlab源程序（见附录）

(2)Matlab仿真图形（三种扰动下的响应曲线）

A 阶跃信号下的响应曲线

B 方波信号下的响应曲线

C 正弦信号下的响应曲线

（3）实际系统的响应曲线

当Q1=500，Q2=700

A 锯齿波信号下的实际响应曲线

B 方波信号下的实际响应曲线

当Q1=300，Q2=500

A锯齿波信号下的实际响应曲线