计算机仿真大赛一等奖作品.docx

计算机仿真大赛一等奖作品.docx

《计算机仿真大赛一等奖作品.docx》由会员分享，可在线阅读，更多相关《计算机仿真大赛一等奖作品.docx（20页珍藏版）》请在冰豆网上搜索。

计算机仿真大赛一等奖作品

编号：

重庆邮电大学

第五届计算机仿真大赛一等奖作品

题号:

C组别：

高年级组

2011年4月6日

单级移动式倒立摆动态控制仿真

摘要

单级移动式倒立摆系统是一个典型的单输入双输出的自然不稳定系统，具有非线性、强耦合、多变量、高阶次等特性，作为控制系统的被控对象，它是一个理想的教学实验设备，许多抽象的控制概念都可以通过倒立摆直观地表现出来。

应用上，倒立摆广泛应用于控制理论研究、航空航天控制、机器人、杂技项杆表演等领域，在自动化领域中具有重要的价值。

研究倒立摆的精确控制对工业生产中复杂对象的控制有着重要的应用价值，因此倒立摆仿真或实物控制实验是控制领域中用来检验某种控制理论或方法的典型方案。

许多抽象的控制概念如控制系统的稳定性、可控性、可观性、系统收敛速度和系统抗干扰能力等，都可以通过倒立摆系统直观的表现出来。

在控制过程中，它能有效地反映诸如可镇定性、鲁棒性、随动性以及跟踪等许多控制中的关键问题，是检验各种控制理论的理想模型。

因而其研究具有重大的理论和实践意义。

本文以单级移动式倒立摆为被控对象，参考古典控制理论设计控制器（如PID控制器）的设计方法和现代控制理论设计控制器（最优控制）的设计方法。

针对单级移动式倒立摆系统的单输入双输出、非线性、强耦合的不稳定性，我们采用线性二次最优LQR控制和PID控制，通过MATLAB软件进行实时控制（摆杆的角度和小车的位移）和系统仿真实验实现对小车位置和摆杆偏角的同时闭环控制，达到了预期效果（动态仿真）。

关键词：

非线性倒立摆实时控制系统仿真动态

1.问题重述

单级移动式倒立摆原理如下图1-1所示，

图1-1

其中小车质量M=3KG;摆杆质量m=0.1KG;摆杆为均匀的，质心l=0.5米;重力加速度g=9.81;设x1为摆杆角度，x2为摆杆角速度，x3为小车位移，x4为小车速度，其中小车水平受力为输入量，建立其状态方程：

采用状态反馈进行控制，其控制原理如下图1-2所示：

图1-2

计算得k1=-446.7287,k2=-91.0448,k3=-244.6376,k4=-110.0896时系统稳定。

建立倒立摆控制仿真系统，当输入信号为阶跃信号和随机信号时，实时演示倒立摆的工作情况。

2.单级倒立摆模型的建立

2.1倒立摆概述单级倒立摆是一个自然不稳定系统，是处于倒置不稳定状态，人为控制使其处于

动态平衡的一种摆，一般是由一个可以在水平轨道上自由移动的小车和倒置摆铰链而成，用一种强有力的控制方法对小车的速度作适当的控制，从而使全部摆杆倒置稳定于正上方，这个系统也叫自动平衡车。

悬挂式倒立摆系统开始工作时，首先使小车按摆杆的自由振荡频率摆动，摆杆随之大幅度摆动，当摆杆接近于倒立摆垂直位置时，自动转换控制方法，使其稳定于倒置状态。

随着倒立摆系统的控制研究的不断深入，倒立摆系统的种类也由简单的单线倒立摆发展为多种形式的倒立摆。

常见的倒立摆系统一般由小车和摆杆两部分构成（另有旋转式倒立摆等形式），其中摆杆可能是一级、两级甚至多级（级数越多，控制难度越大），其长度也可能是可变化的。

控制的目标一般都是通过给小车施加一个水平方向的力，使小车在期望位置上稳定，而摆杆达到竖直向上的动态平衡状态。

2.2数学准备对系统建立数学模型是系统分析、设计的前提，而一个准确又简练的数学模型将

大大简化后期的工作。

为了简化系统分析，在实际的模型建立过程中，要忽略空气流

动阻力，以及各种次要的摩擦阻力。

这样，可将倒立摆系统抽象成小车和匀质刚性杆组成的系统，如下图所示。

本系统内部各相关参数定义如下：

M小车质量

m摆杆质量

b小车摩擦系数

l摆杆转动轴心到杆质心的长度

I摆杆惯量

f加在小车上的力

x小车位置

φ摆杆与垂直向上方向的夹角

θ摆杆与垂直向下方向的夹角（考虑到摆杆初始位置为竖直向下）下图是系统中小车和摆杆的受力分析图。

其中，N和P为小车与摆杆相互作用力

的水平和垂直方向的分量。

注意：

在实际倒立摆系统中检测和执行装置的正负方向已经完全确定，因而矢量方向定义如图，图示方向为矢量正方向。

应用Newton方法来建立系统的动力学方程过程如下：

分析小车水平方向所受的合力，可以得到以下方程：

Mx=

f-bx-N

由摆杆水平方向的受力进行分析可以得到下面等式：

2

N=md

（x+lsinθ）

dt2

即：

N=mx+mlθcosθ-mlθ2sinθ

把这个等式代入上式中，就得到系统的第一个运动方程：

（M+m）x+bx+mlθcosθ-mlθ2sinθ=F

（2-1）

为了推出系统的第二个运动方程，我们对摆杆垂直方向上的合力进行分析，可以得到下面方程：

2

P-mg=-md

（lcosθ）

dt2

即：

P-mg=mlθsinθ+mlθ2cosθ

力矩平衡方程如下：

-Plsinθ-Nlcosθ=Iθ

注意：

此方程中力矩的方向，由于θ=π+φ,cosφ=-cosθ,sinφ=-sinθ，故等式

前面有负号。

合并这两个方程，约去P和N，得到第二个运动方程：

（I+ml2）θ+mglsinθ=-mlxcosθ

（2-2）

2.3模型建立

2.3.1微分方程建立

设θ=π+φ，当摆杆与垂直向上方向之间的夹角φ与1（单位是弧度）相比很小，

即φ<<1

时，则可以进行近似处理：

cosθ=-1，sinθ=-φ，（dθ）2=0。

为了与控

dt

制理论的表达习惯相统一，即u一般表示控制量，用u来代表被控对象的输入力F，线性化后得到该系统数学模型的微分方程表达式：

⎨

⎪⎧（I+ml2）φ-mglφ=mlx

⎪⎩（M+m）x+bx-mlφ=u

（2-3）

2.3.2传递函数建立

对方程组（2-3）进行拉普拉斯变换，得到

⎪⎧（I+ml2）Φ（s）s2-mglΦ（s）=mlX（s）s2

⎨

⎪⎩（M+m）X（s）s2+bX（s）s-mlΦ（s）s2=U（s）

注意：

推导传递函数时假设初始条件为0。

由于输出为角度φ，求解方程组（2-4）的第一个方程，可以得到

（2-4）

2

X（s）=[（I+ml）-

ml

g]Φ（s）

s2

把上式代入方程组（2-4）的第二个方程，得到

⎡+2⎤

⎡2⎤

（I

（M+m）⎢

ml）-

g（I+ml

⎥Φ（s）s2+b⎢

）g

-⎥Φ（s）s-mlΦ（s）s2=U（s）

⎣ml

s2⎦

⎣ml

s2⎦

整理后得到以输入力u为输入量，以摆杆摆角φ为输出量的传递函数：

mls2

2

G（s）=Φ（s）=q

U（s）

s4+b（I+ml）s3-（M+m）mgls2-bmgls

qqq

其中q=[（M+m）（I+ml2）-（ml）2]

若取小车位移为输出量，可得传递函数：

（I+ml2）mgls2-

2

2

G（s）=X（s）=qq

U（s）

s4+b（I+ml）s3-（M+m）mgls2-bmgls

qqq

2.4状态空间模型

2.4.1状态空间模型的建立状态空间模型的一般形式如下：

x=

Ax+Bf

y=Cx+Df

2.4.2系统向量的设定

设x1为摆杆角度，x2为摆杆角速度，x3为小车位移，x4为小车速度

T

x=

⎡⎣x1

x2

x3

x4⎤⎦

0

0

1

0

0

1

0

0

0

T

0]

⎡10⎤

⎢00⎥

⎡

0

1

0

0⎤

⎢M+m

⎢

g0

0

⎥

0⎥

A=⎢Ml

⎥C=⎢⎥

⎢0001⎥

⎢00⎥

⎢⎥

⎢-mg

000⎥

⎢01⎥

⎣⎢M⎥⎦⎣⎦

B=⎡0-1

01⎤

D=[00

T

⎣⎢MlM⎥⎦

2.4.3系统的状态空间模型

将上面的向量带入状态方程的一般形式，得到如下的状态方程：

3.单级倒立摆模型性能分析

3.1系统能控性分析

系统能控性就是研究系统的全部状态是否都会受到输入的影响，从而实现对系统状态的控制。

能控性体现了系统输入对状态的控制能力。

[1]

根据系统能控性定义可知，系统完全能控的充分必要条件是：

能控性矩阵

Wc=⎡⎣AAB

An-1B⎤⎦为行满秩。

仿真程序：

[2]n=length（A）;Q=zeros（n）;Q（:

1）=B;

fori=2:

n%按照公式Q=[BABA^2*B....A^（n-1）*B]Q（:

i）=A*Q（:

i-1）;

endm=rank（Q）;if（m==n）

disp（'该系统可控！

'）;

melse

disp（'该系统不可控！

'）;

mend

仿真结果：

仿真结果分析：

图3-1-1

由仿真结果知，可控性矩阵的秩为4，即为行满秩，故该系统具有可控性。

3.2系统能观性分析系统的能观性就是系统状态变量的任何运动完全可以由输出来反映。

能观性表征了

输出对状态的反应能力。

[1]

根据系统能观性定义知，能观性完全能观的充分必要条件是：

能观性矩阵

W0=⎡⎣CCA

T

CAn-1⎤⎦

为列满秩。

仿真程序：

V=obsv（A,C）%求能观性矩阵rank（V）%求能观性矩阵的秩

仿真结果:

图3-2-1

由仿真结果知，可观性矩阵的秩为4，即为为列满秩，故该系统具有可观性。

3.3阶跃响应分析

3.3.1开环系统的阶跃响应仿真程序：

%求解系统的特征值%p=eig（A）t1=0:

0.01:

1;

step（A,B,C,D,1,t1）%阶跃响应

仿真结果：

图3-3-1原系统阶跃响应曲线结果分析：

由于特征方程有一个正实部的根，故系统不稳定。

3.3.2闭环系统的阶跃响应

加入反馈K=[-446.7287-91.0448-244.6376-110.0896]后的阶跃响应

仿真程序：

%带反馈T1=0:

0.01:

5;u1=0.25*ones（size（T1））;[Y,X]=lsim（A-B*K,B,C,D,u1,T1）;plot（T1,Y（:

2））;

title（'加反馈后系统的阶跃响应曲线'）;

gridon

仿真结果：

仿真结果分析

图3-2-2

由图3-2-2可知，系统的快速性很好，过渡过程时间不超过2.5秒。

相比开环系统的阶跃响应（图3-3-1）来说，闭环系统的稳定性能明显较好。

4.单级倒立摆模型动态仿真

4.1小车模型的建立

4.1.1车身的创建

如图4-1-1创建的简易车身

图4-1-1

4.1.2车轮的创建

根据圆的绘图公式x=a+r*sin（t）;y=b+r*cos（t）;

如图4-1-2所示创建的车轮

figure

（1）

图4-1-2

plot（x5-4,y5,'ko',x5+4,y5,'ko',[x1-4,x2+4],[y1-1,y1-1],'k-',[x1-4,x2+4]

[y2+1,y2+1],'k-',[x1-4,x1-4],[y1-1,y2+1],'k-',[x2+4,x2+4],[y1-1,y2+1],

'k-'）

axis（[-1515-510]）

4.4.3整车模型的创建

4.2倒立摆模型的建立

4.2.1摆杆和铰链的创建

4.3单级倒立摆动态仿真

4.3.1输入信号为阶跃信号的动态显示阶跃信号的产生[4]

%阶跃函数

t0=0;t1=-2;t2=10;

t=t1:

0.001:

t2;

x=（t-t0）>=0;

plot（t,x,'r-'）,axis（[t1t2-0.52]）;xlabel（'t'）;ylabel（'x'）;title（'阶跃信号'）;gridon

4.3.2输入信号随机信号时的动态显示

%随机信号的产生[5]Fs=500;%频率

x=-1+（1-（-1））.*rand（100,1）;%从-1到1创建随机信号

plot（x）;

title（'随机信号'）;

图4-3-2

4.3.3系统对阶跃信号和随机信号的动态响应

图4-3-1图4-3-2

动态仿真描述：

图4-3-3

开始给单级移动式倒立摆一个阶跃信号（随机信号），倒立摆随之摆动，倒立摆的角度和小车位移就会发生变化，通过反馈控制，给小车施加外力使小车产生一定的加速度，从而使倒立摆基本上趋于稳定。

5.结语

计算机仿真是提高学生综合素质的一种有效方式，在这次竞赛过程中，我们分工协作，相互配合，相互鼓励，不但学会了matlab仿真的相关知识，了解到一些关于控制理论的算法，而且实现了单级移动式倒立摆的动态控制，理论联系实际，让我们受益非浅。

附录

%【小车动态仿真源程序】smallcarmove.m

%设x1为摆杆角度，x2为摆杆角速度，x3为小车位移，x4为小车速度M=3;%小车质量

m=0.1;%摆杆质量l=0.5;%质心g=9.81;%重力加速度k=5;

k1=-446.7287;%x1反馈系数k2=-91.0448;%x2反馈系数k3=-244.6376;%x3反馈系数k4=-110.0896;%x4反馈系数

%状态方程的建立a21=（M+m）*g/（M*l）;a41=-m*g/M;

A=[0100;a21000;0001;a41000];B21=-1/（M*l）;

B41=1/M;B=[0;B21;0;B41];

C=[1000;0100;0010;0001];%输出为小车位移和摆杆角度D=[0;0;0;0];

K=[k1k2k3k4];

dt=0.01;

%系统初始值zf=[0010]';z=[0020]';figure

（2）

fori=1:

500

z1=（A-B*K）*（z-zf）;z=z+z1*dt;y（:

i）=z;

%作动画

x1=z（3）-k/2+z（4）*dt*i;x2=z（3）+k/2+z（4）*dt*i;y1=4;y2=0;r1=2;x4=（x1+x2）/2+r1*sin（z

（1））;y4=y2+r1*cos（z

（1））;t=0:

pi/100:

2*pi;x3=x4+r1*cos（t）;y3=y4+r1*sin（t）;r2=0.5;%车轮半径x5=x1+r2*cos（t）+2.5;y5=y1+r2*sin（t）-3;r3=0.1;%铰链半径x6=2.5+r3*cos（t）;y6=3+r3*sin（t）;

%%%%车前轮车后轮铰链杆地平线上边下边左边右边%%%

plot（x5-4,y5,'ko',x5+4,y5,'ko',x5,y6,'c-',[（x1+x2）/2,x4],[y2+3,y4+3],'r-',[-20,20],[0.4,0.4],[x1-4,x2+4],[y1-1,y1-1],'k-',[x1-4,x2+4],[y2+1,y2+1],'k

-',[x1-4,x1-4],[y1-1,y2+1],'k-',[x2+4,x2+4],[y1-1,y2+1],'k-'）

axis（[-1515-510]）axissquaretitle（'单级移动式倒立摆动态模型'）pause（0.001）;

if（i==1）

pause（0.01）;

end

if（abs（z

（1）-zf

（1））==0）&（abs（z（3）-zf（3））==0）%判定程序结束

break;

endend

参考文献

[1]王平，谢昊飞，蒋建春，计算机控制技术，北京：

机械工业出版社，2010。

[2]冯辉宗，岑明，张开碧，控制系统仿真，北京：

人民邮电出版社，2009。

[3]魏巍，控制工程工具箱技术手册，北京：

国防工业出版社，2004。

[4]刘正君，MATLAB科学计算与可视化仿真宝典，北京：

电子工业出版社，2009。

[5]王丽，信号产生，

2011年4月5日。