基于matlab在连续性系统中的仿真研究.docx
《基于matlab在连续性系统中的仿真研究.docx》由会员分享，可在线阅读，更多相关《基于matlab在连续性系统中的仿真研究.docx（22页珍藏版）》请在冰豆网上搜索。
基于matlab在连续性系统中的仿真研究
基于matlab在连续性系统中的仿真研究
摘要：
从数学建模到仿真系统的搭建，再到加进控制环节进行实时控制，最后得出结果的过程中，参考了大量的资料，通过对比整合，设计出了适合自己的一套实验方法：
本研究的核心在于构建相应的倒立摆数学模型，首席按通过牛顿•欧拉法完成相关模型的基本搭建，之后通过动态系统空间状态方程来确定其系数矩阵的值，完成上述匸作之后，通过MATLAB编程的方式，将构建的数学模型转化成相应的传递函数，从而获得传递函数模型。
此后便可在此基础之上，通过Siinuluik软件通过PLD控制算法实现相应的仿真工作，先用连续系统的设计方法设讣出模拟控制器，然后在满足一定条件下，对其进行离散化处理，（采用加零阶保持器的Z变换法）形成数字控制器。
接着进行PID参数整定，利用试凑法，根据PED控制器各部分组件和性能之间的联系，在某个初始化参数的条件下重复进行试凑工作，最终得到较好的结果。
研究过程中，系统的控制非常稳定，性能较好。
关键词：
Matlab；仿真；PID；连续系统
1绪论1
1・1连续性系统简介1
1.2连续性系统基本处理方法1
1.2.1数值积分法1
1.2.2替换法2
1.2.3离散相似法3
1.3MATLAB在连续性系统的应用5
2MATLAB基础6
2.1MATLAB语言介绍6
2.2MATLAB工具箱6
2.3MATLAB对连续性系统的仿真性能分析错误!
未定义书签。
3直线一级倒立摆系统设计7
3.1倒立摆系统概述7
3.2倒立摆系统工作原理7
3.3倒立摆系统数学模型建立8
3.3.1确定系统输入输出量及中间变量：
9
332受力分析，列写运动方程：
9
333系统非线性方程的线性化：
10
3.3.4零初始条件下的拉氏变换：
10
3.3.5代入参数求解传递函数：
10
4基于MATLAB的一级倒立摆系统仿真结果及分析11
4.1校正前系统性能仿真分析11
4.1.1稳定性分析：
11
4.1.2阶跃响应分析：
11
4.1.3频率特性分析：
12
4.2PID控制设计校正装置13
421计算开环增益K：
13
4.2.2计算校正后的增益交界频率3C：
15
423计算参数确定传递函数：
15
4.3控制实验仿真分析15
5总结和展望20
参考文献21
1绪论
1.1连续性系统简介
所谓连续系统，其实就是体系的状态会在时间迁移的条件下产生光滑且连续的改变的一个系统。
其也包含那些因离散采样而引起非连续改变的体系。
这类系统通常能够通过微分方程组进行表示。
如果方程之中的相关系数皆是常数，那么可将之称作定常系统，如果相关系数会因时间的改变发生改变，那么可将之称作。
除此之外，相关的数学模型也可分为多个类别，除了连续模型（即通过微分方程进行表示）之外，还有离散模型（通过差分方程进行表示）以及二者的混合模型。
如果某个体系能够通过非线性的常微分方程组进行表示，那么其就属于非线性系统。
其中的各个方程不仅含有状态变量及其各阶导数的线性组合，也含有非线性项。
其非线性特性主要表现在：
（1）线性系统同时满足均匀性与叠加性。
线性方程必定有解，并且其解为零输人响应和零状态响应之和。
非线性系统则不具有这些特征，即对非线性系统，通常叠加性原理不成立，非线性方程通常没有解析表达形式的解。
（2）给定一个线性系统，对于不同的初始状态，它只有一种类型的运动。
但对于非线性定常自山系统则不然，初始状态不同，它对应的运动类型可以不一样。
非线性系统有三种类型的运动，即稳定运动、不稳定运动和周期运动。
1.2连续性系统基本处理方法
1.2.1数值积分法
如果通过相关的状态方程来表示被模拟的体系，其对应的一阶微分方程以及初始值可以通过下式表示：
（1-1）
yM=yo
在这之中，Y代表维度为n的状态向量，而F则代表维度为n的向量函数。
对于上式，当r=…,U+i时其对应的连续解可以写成：
ln*l'（H-1
Y（tn+J=Y（t°）+JF（⑴〃=g）+jF（f,y）df（1.2）
*0l0
设rn=r（^）,令
却=監+0（1-3）
则有：
匕+i=Y（J+J
也就是说，
a,«jF（t,Y）dt（1-4）
Zn
如果y”准确解”/”）为近似值，g是准确积分值的近似值，则式（1-4）就是式（1-2）的近似公式。
换句话说，连续系统的数值解就转化为求两个相邻时刻之间的积分。
综上所述，（1・1）的数值求解策略本质上就是求得该式的实际解在相应的离散时间点人“2・・叫〈…之中的近似值E，E，…必…，对于处在相邻时刻的两个点，其间距hn=t^-tn,通常被命名为步长，一般令h”=〃，保持不变。
1.2.2替换法
以上使用的通过数值积分对连续体系进行模拟的策略尽管理论上较为完善、求解精度也较为不错，然而过程较为繁琐，需要进行大量的讣算，一般仅将之应用于对于效率没有过高要求的模拟工作之中。
如果对于模拟工作有着实时性的需要，要求计算速度快，以便能在一个采样周期内完成全部计算任务，这就必须使用其他的高效率的求解方法。
通过数值积分策略实现对连续体系的模拟，本质上已对其实施了离散化的操作，只是在进行离散化的各个步骤之中都会使用连续体系的模型，每一步都是离散之后再进行求解。
因此，可以先将相关连续模型实施离散化操作，获得等价的离散模型，在这之后的各个计算工作中可以不必再涉及最初的连续模型，而是能够直接使用该离散模型。
通过对连续模型简化得到相应的离散模型，能够很好地适应计算机上的计算环境，不但能够应用在连续体系的数字模拟工作中，同时也能够通过计算机实现相应的数字控制器。
替换法的核心思路为：
对某个传递函数G（s）,寻找将s域映射至z域的策略，从而实现属于S域的变量到z平面之间的映射工作，从而能够将和连续体系相关的G（s）转化为离散的G（z）o之后可以通过G（z）进行相应的z反变换，得到对应的通过差分方程进行表示的离散模型，从而实现高效率的计算工作。
根据z变换理论，s域到z域的最基本的映射关系是或s=丄In?
T
如果按这一映射关系直接代入G（s）,得到的G（z）是相当复杂的，不便于算法实现，所以往往借助于Z变换的基本映射关系Z=0或$作一些
T
简化和近似处理。
1.2.3离散相似法
离散相似法能够把某个连续的体系实施离散化操作，并获得与之对应的差分方程形式的离散模型。
得到对应的离散模型的策略包括两类：
（1）对传递函数实施相应的离散化操作，从而获得离散化的传递函数一—这类模型通常叫做“频域离散相似模型”；
（2）对状态方程实施相应的离散化操作——这类模型通常叫做“时域离散相似模型⑶”；
为了实现连续体系的数字模拟工作，首先可向体系之中添加相应的虚拟采样器以及保持器，具体的架构可以参考下图：
图1-1连续系统离散化结构图
附注：
上图之中给出的采样开关以及保持器实际上是不存在的，而是为了将（1・5）式离散化而虚构的。
然后通过z变换实现对应体系脉冲传递函数的求解，之后可在此基础之上获得体系所对应的差分方程。
根据上图不难得到相应的脉冲传递函数模型，其可通过下式进行表示：
在这之中，G力G）代表保持器的传递函数，并且可根据保持器的种类获得相应的G⑵，具体可参考下表：
表1・1不同保持器的G⑵
保持器的传递函数Gh（s）脉冲传递函数G⑵
x=Ax+Bu（1-6）
如果在体系的输入以及输岀部分都添加相应的采样控制开关，在输入部分额外添加保持器以确保相关信号能够良好地还原，那么此时可获得下图所给岀的结构：
图1・2采样控制系统结构图
若对方程（1・6）式两边进行拉普拉斯变换，得：
sX（s）-X（0）=AX（s）+BU（s）
即：
（N-A）X⑶=X（0）+BU（s）
以（si-A）'1左乘上式的两边可得：
X（$）=（si-A）-1X（0）+（$/（1-7）
考虑到状态转移矩阵：
==rl（1-8）
故对（1・7）式反变换可得：
X（0=eAtx（O）+矗从~r）BU（r）dT（1-9）
此为（1・5）式的连续解，山此可推导出系统的离散解。
根据上式，n及ii+l两个相连的采样瞬间，有：
X（kT）=eAkTx（Q）+\^eA（X[伙+1）T]=eA（k+叭（0）+瞪+”尹伙+V）T~^BU（T）dr（1-11）
将（1・10）式减去（1・10）式后乘以評，得：
X[伙+1）T]=eATX（kT）+球+1）7/[伙+l）T~r^BU（T）dT（1"2）
将（1-11）式右边积分进行变量代换，即令：
r=kT+t（1-13）
则得：
X[（k+1）7]=eATX（kT）+eA（T~{）BU（kT+t）dt（1-14）
但由图1・2可知：
若系统采用零阶保持器时，则两个采样点之间输入量可看做常数，即u（nT+t）=u（nT）,这样（1-14）式可写为:
X[伙+1）T]=eATX（kT）+^eA（T~l）BdtU（kT）
=G（T）X（kT）+H（T）U（kT）
G（T）=eA1
式中•H（T）=^eA（I
1.3MATLAB在连续性系统的应用
瑞典的著名高校Limds大学曾经开展过时间为3个月的极地研究工作，并在这个过程中大量运用到了MATLAB和相关工具箱，其主要用于考察北极的放射性物质和地环境之间的作用关系。
相关研究者通过MATLAB来完成数据处理以及分析匸作，同时通过提供的神经网络工具箱实现了对北冰洋之上漂浮木材的年轮图像识别工作，从而确定相关浮木的来源。
除此之外，Forsmark核电站也通aMATLAB来完成对核反应堆的功率优化工作，相关的研究者能够在各个堆芯之中得到海量的数据，并在此基础上求解各燃料棒以及控制设备坐标的最优值,从而获得最高的功率。
这个数据处理过程的规模非常大，且复朵度也极舟，涉及到超过1700个结点的处理工作，为能精简这个流程，Forsmark通过MATLAB开发得到了相应的GUI交互界面，其能够让毫无经验的人员进行数据的处理和分析工作。
不仅如此，Forsmark还进一步通过MATLAB对一系列失效和扰动的场景进行建模处理。
一旦反应堆出现了扰动现象，应当对相关信息进行详实的分析，从而实现引发扰动因素的定位工作，利用MATLAB的种种优势，上述工作所耗费的时间已经从原本的7天降低到15分钟。
Calspan先进技术中心曾对相关的飞行器进行仿真测试，并有效提高了其性能。
相关的模拟软件工作与一系列并行的浮点DSP中，硬件部分则使用的dSPACEo在地面阶段，首先通过Simulink模块实现相关飞行器及其控制系统的建模以及模拟工作，之后通过Workshop产生对应的C语言代码，并传输到飞行器的DSP之中。
进行试飞之时，试飞员能够通过DSP模块对该系统进行考察,E行过程之中涉及到的一系列参数也能够实时调整，同时也可下载到硬件之中进行测试，从而有效降低了成本，提高了效率，能够在控制系统原型的实际构建之前便将相关试验全部完成。
2MATLAB基础
2.1MATLAB语言介绍
MATLAB诞生初期，其仅仅是用于为数据分析工具Eispack与Linpack提供接口，其最早通过Fortran实现，之后通过C实现。
1975年左右，M公司将其推出面向市场，并且将其功能逐渐丰富，开发了文字图像的处理、数值的运算以及符号的解译等功能，将其编程语言升级为面向对象语言，并以此为用户界面的设计语言，从而将其打造成了一个能够在多个领域、多个学科都很好地发挥功能的软件。
其版本也逐渐的更新换代，其有效综合了各种数据处理功能，包括矩阵求解、数据可视化、非线性体系的分析等等，SJ询在工程领域以及科研工作之中都得到了普遍的应用，其主要包括下面儿个部分：
2.2MATLAB工具箱
现阶段,MATLAB已经拥有了超过37个工具箱，同时其主要包括两个类别，分别是功能类以及领域类。
后者通常依赖于较高的特定领域专业知识，主要有信号处理、自控、逻辑分析、图象处理、通讯、偏微分方程等等很多细化的工具箱，可应用的领域涵盖了化学、数学、工程还有经济等等很多的领域范圉。
功能型的工具箱主要是增强其图形文字处理能力、模型建立和仿真功能、数值符号运算解议以及和软硬件系统的通信等一系列作用。
使用者可以直接对相关工具箱的源代码进行自行定制，以满足特定的需求。
3直线一级倒立系统设计
3.1倒立摆系统概述
倒立摆系统是用于进行控制理论分析工作的一个非常经典的设备，其有着非常精简的构造，十分低廉的成本，同时相关参数以及构造可以灵活地更改。
该系统处于非稳态，涉及变量多、阶次高、耦合高。
通过该系统，可以有效体现出控制流程里面的可能问题，包括稳健性、可定位性等等，该系统能够用于对相关的控制理论进行检验。
在各类控制理论的基础之上，相关研究者已经设计了各种类型的倒立摆系统。
3.2倒立摆系统工作原理
图3-1便携式一级倒立摆试验系统总体结构图
该体系总的来说山摆杆、支架、导轨、驱动设备以及小车等部分构成。
主体、驱动器、电源和数据采集卡都置于实验箱内，实验箱通过一条USB数据线与上位机进行数据交换，另有一条线接220v交流电源。
便携式直线一级倒立摆的工作原理如图3-2所示：
据集数采卡
编码器
编码器*
图3-2便携式一级倒立摆工作原理图
数据采集卡能够获得各个编码器产生的数据，其中，旋转编码器和摆杆之间处在相同的轴线之上，而电动机和车辆之间以皮带相连，因此能够求出摆杆的偏角和车辆的偏移，并可进一步得到角速度以及线速度，之后就能够利用相关的控制理论求得对应的控制量。
后者会山讣算机传输到驱动设备之中，从而能够进行电机的控制工作，并保证摆杆处在竖直状态。
该闭环体系拥有反馈功能，主要的控制U的是在静止和运动状态下实现稳摆。
如果输入0g=O,则此时的偏差=0,控制设备并不会执行操作；否则，偏差不等
于0,控制设备会进行分析处理，并对电机进行控制，引起摆杆的运动，获得的摆角通过角位移传感器作用于输出量，达到减小偏差的H的。
根据控制原理绘制出控制方框图如图3・3所示：
角位移传感器
图3-3便携式一级倒立摆控制原理方框图
3.3倒立摆系统数学模型建立
该体系主要包括小车以及摆杆两个模块。
摆杆能够在垂直的平面内自山地摆动，小车能够在轨道上自山移动。
忽略摩擦力以及空气阻力，不妨将该系统视作均质杆体和小车构成的一个刚体，因此可在经典力学的基础之上得到相应的动力模型，对其进行受力分析可以获得对应的数学模型，在平面坐标系中该系统的结构如图4・1所示，规定逆时针方向的转角和力矩均为正。
模型参数符号如下：
M:
小车质量；m:
摆杆质量；/:
摆杆转动惯量；I:
摆杆质心到转轴的距离;
L:
摆杆长度；F:
外加在小车上的力；N:
摆杆给小车的力；P:
小车给摆杆的支持力；0：
转角；%：
小车位置；b:
小车与导轨间的摩擦系数。
3.3.1确定系统输入输出量及中间变量:
输入量：
F:
加在小车上的力
输出量：
0:
转角
中间变量：
尢：
小车位置
3・3・2受力分析，列写运动方程：
在水平方向上，根据上图可以得到小车的受力为：
F-bx-N=Mx
（1）
在水平方向上，根据上图可以得到摆杆的受力为：
N=m—-（x一Isin0）
（2）
dr
2
根据式
（2）有：
N=ma-ml（/）cos（j）+ml（j>sin0（3）
综合以上各式，能够获得体系的运动方程1:
2
+in）x+bx-ml（j）cos（j）+ml（psin°=F（4）
在竖直方向上，根据上图可以得到摆杆的受力为：
d?
P一mg=m（/cos0）（5）
由（5〉得：
P-mg=-mlsin-ml（/）cos0（6）
对摆杆进行力矩分析，有：
I（/）=Pls\n（/>+Nlcos°（7）
联立式（3）、（6）以及（7）,可以获得体系的运动方程2：
（/+m卩）（/>一mglsin（/>=mlxcos0（8）
因此可以获得体系的2个运动方程，如下所示：
（2
I（M+m）x+hx-ml（/）cos（f>+ml（f>sin（/）=F（（/+ml2）。
一mglsin°=mlxcos（j>
3・3・3系统非线性方程的线性化：
对于该系统，摆杆通常不会产生较高的偏角，也就是说6«1,因此能够近似化处置。
对于该体系，当（甩0°）为（0,0）之时处在平衡态，因此：
sin0=0.cos&=l,=0・
因此可得线性化的方程：
{
（M+m）x+bx-mlQ=F
（/+ml2）0-=mlx
33.4零初始条件下的拉氏变换：
I
（M4-（s）+bsX（5）-mis20（f）=F（s）
（/+ml2）s2^（s）-=mlX（s）s2
3・3・5代入参数求解传递函数:
该体系的物理参数为:
小车的质量M：
0.618kg
摆杆的长度L：
0.350m
摆杆的质量th：
0.0737kg
摆杆质心到转轴距离h0.1225m
把相关数据回代到上式之中，即可计算岀摆杆的偏角和车辆偏移之间的传递
函数:
0（s）6.122s2
X（s）_s-60
4基于MATLAB的一级倒立摆系统仿真结果及分析
4.1校正前系统性能仿真分析
4・1・1稳定性分析：
可通过下面的代码得到体系进行矫正之询的Nyquist图，具体结果可参考下图：
num=[6.122];
den=[10-60];
nyquist（num,den）
图4-1校正前系统开环Nyquist图
s_xvAjeu一6PUJ_
由图分析知，系统开环右极点数P=ItN=0。
此时Y/Z=P+JV=1,所以体系处在非稳态。
4.1.2阶跃响应分析：
可通过下面的代码得到体系进行矫正之询的阶跃响应图，具体结果可参考下图：
G=6.122/（sA2-60）;
GB=G（1+G）;step（GB）
.5
1
15
6
0>pm=dE<
23456
Time（seconds）
图4・2校正前系统阶跃响应曲线图
根据该结果不难发现，体系的阶跃响应数据不会随时间迁移而降低，并处在逐渐发散的状态中，所以体系处在非稳态。
4.1.3频率特性分析：
可通过下面的代码得到体系进行矫正之前的Bode图，具体结果可参考下图:
num=[6.122];
den=[l0-60];
bode（num,den）;
BodeDiagram
ooO
-2-4-6
mp）apmc6ew
-80
-179
-179.5
・180
-180.5
-181
10°101102
Frequency（rad/s）
图4・3校正前系统Bode图
山图分析知，相频特性图中相位裕度为0,因此体系处在非稳态，可以通过频率控制器进行调整。
4.2PID控制设计校正装置
构建相应的控制器，从而让校正之后体系的稳态位置误差系数Kp=10,同
时对应的相位裕度Y（Gc）=50°o
在进行校正工作之前，体系的传递函数可通过下式表示:
如果用K代表开环增益，那么进行校正之后对应的传递函数可通过下式表示：
421计算开环增益K：
由于匕=10■因此!
乌G&（s）=10：
，、八Ts+l6.122
KaGc（s）G（s）=K——pr—777=
“+产一60
a
可以求岀:
K=98
因此:
可通过下面的代码得到KG（s）的Bode图，具体结果可参考下图:
num=[6.122*98];
den=[l0-60）;
bode（num,den）
图4-4添加增益后的便携式倒立摆的Bode图
不难发现，此时体系的相位裕量依旧是0°,因此对于相关的矫正设备，其应出产生的最大相位超前角监可通过下式进行讣算
0加=人一升+[0（呎）一0（"；）1
对本系统:
血=50>—0+[-180°-（一1）]=50>
通过0m可以得到相应的⑴
4.2.2计算校正后的增益交界频率
为能充分显示出相关校正设备的相位超前效果，可以把校正设备的最大超前相角确定为体系增益交界频率也就是3尬=Sc，在该处，冇：
9（%）+厶（3』=0
带入可得：
201g\/a4-厶（3』=0
可以计算岀:
3c=39.Sradfs
4.2.3计算参数确定传递函数：
y[a
com=o）c==39.5rad/s
得：
T=0.06965
则校正后系统的开环传递函数为：
、0・0696s+16.122
K叫（s）G（s）=985092S+2-60
4.3控制实验仿真分析
首先确定临界增益心和临界振荡周期几，只加入P控制器，令Kp=9得:
可通过下面的代码得到体系进行增益之后的阶跃响应图，具体结果可参考下图：
s=tf（U）
G=9*6.122/（sA2-60）;
GB=G（1+G）;step（GB,10）
10
012345678910
Time（seends）
图4・5Kp=9时系统的阶跃响应曲线图
根据该结果不难发现，体系的阶跃响应数据不会随时间迁移而降低，并处在逐渐发散的状态中，所以体系处在非稳态。
因此提高P控制量，令Kp=100得系统的阶跃响应曲线图如图4-6：
StepResponse
图4・6Kp=100时系统的阶跃响应曲线图
根据该结果不难发现，体系的阶跃响应数据处在逐步收敛的状态中，为让该曲线处在震荡状态，应降低Kp,故令心=90得系统的阶跃响应曲线图如图4-7：
25
StepResponse
图4・7知=90时系统的阶跃响应曲线图
根据该结果不难发现，临界的增益K“=90,此时的周期几=0.4。
因此可算岀Kp.八以及丁小
Kp=0・6Ku=0.6X90=54
7\=0.5Tu=0.2
Td=0.125几=0.05
z、16.122
G*（s）=54x（1+0.05sH）x—
0.2ss2-60
可通过下面的代码得到体系进行PED校正之后的阶跃响应图，具体结果可参考下图：
s=tf（U）
G=54*（l+0.0375*s+l/0.15/s）*6.122/（sA2-60）;
GB=G（1+G）;step（GB,4）
2
1.8
1.6
1.4
a1.2
0.8
0.6
0.4
0.2
005115225335
Time（seconds）
图4-8心=90时PID校正后系统的阶跃响应曲线图
根据该结果不难发现，曲线处在剧烈震荡之中，所以要提高比例系数，取心=100可绘制出PLD校正后系统阶跃响应曲线
升级会员 