基于Simulink状态空间建模的系统分析方法程序实现.docx
《基于Simulink状态空间建模的系统分析方法程序实现.docx》由会员分享，可在线阅读，更多相关《基于Simulink状态空间建模的系统分析方法程序实现.docx（23页珍藏版）》请在冰豆网上搜索。
基于Simulink状态空间建模的系统分析方法程序实现
基于Simulink状态空间建模的系统分析方法程序实现
荆晓莉
（陕西理工学院物电学院电子信息科学与技术1101班，陕西汉中723001）
指导老师：
***
[摘要]无论用何种方法求高阶连续系统解析解都是十分棘手的问题。
实际上，科学研究和工程应用中更多地需要系统的数值解。
调用Matlab的Simulink工具包，可以用模块图标方法来编程，并通过运行系统模型文件的方法直接给出连续系统的数值解，而不需要将连续系统转换为离散系统再求解。
对于复杂LTI系统直接写出微分方程再给出状态空间矩阵很困难，我们获得系统状态空间矩阵的方法是：
先将系统映射到s域，列出解（系统输出）函数的像函数满足的代数方程组，解代数方程组给出系统函数H（s）,再调用Matlab的函数[a,b,c,d]=tf2ss（num,den）就得到系统的状态空间矩阵，从而完成系统的描述，再创建系统模型文件，写入状态空间矩阵、输入信号、初值条件及运行相关参数，最后编程调用sim（）函数运行模型文件给出连续系统的数值解。
[关键字]连续系统;离散化;Simulink;M文件
TheimplementationofthesystemanalysismethodbasedonSimulinkstatespacemodeling
JingXiaoli
（Grade11,Class1,MajorElectronicInformationScienceandTechnologyDepartmentofPhysics,ShannxiUniversityofTechnology,Hanzhong,723001）
Tutor:
LongShuming
AbstractItisadifficultproblemforthehigherordercontinuoussystemtosolvetheproblem.Infact,thenumericalsolutionofthesystemismoreneededforscientificresearchandengineeringapplication.CallingMatlabSimulinktoolkitcanbeprogrammedbyusingthemethodofmoduleicon,andthedirectmethodofoperationsystemmodelfilearegivencontinuoussystemofnumericalsolution,withouttheneedtoconvertthecontinuoussystemtodiscretesystem,tobesolved.Continuousmodel,mappedtodirectlywriteforcomplexLTIsystemdifferentialequationsandgivethestatespacematrixisverydifficult,weobtainthesystemstatespacematrixmethodis:
firstthesystem'sdomain,alistofsolutions（output）functionasfunctionsatisfiesthealgebraicequationgroup,thesolutionofalgebraicequationsgivesthesystemfunctionH（s）,andthencalltheMatlabfunction[A,B,C,D]=tf2ss（num,den）systemstatespacematrix,thuscompletingthesystemdescription,andthencreateasystemfile,writethestatespacematrix,inputsignal,theinitialconditionsandoperatingparameters,programmingatlastcalltosimrunthemodelfileissystemofnumericalsolutions.KeywordsContinuoussystem,Discrete,Simulink,Mfile
最近几年科学不断发展，系统的结构也渐渐复杂。
所以系统变成了一个值得研究的话题，系统输出的变化情况与系统内部一些变量的关系也值得研究。
本课题基于对信号与信息处理课程中用Matlab/Simulink建模及应用分析实际问题的深入研究。
由于Simulink可以对动态系统进行建模，仿真及分析的软件包，依托于Matlab的丰富的仿真资源，可以使用其快速对系统进行仿真，准确分析所描述动态系统[1]。
除此之外，还可用其中的M文件编程以改变系统中参数，更加方便快捷地研究参数对系统的影响，包括求解系统，绘制波形并进行频谱分析等等。
在许多情况下，人们不禁关心系统输出的变化情况，而且还要研究与系统内部一些变量有关的问题，比如，系统的客观性和可控性行。
系统的最优控制与设计等问题。
为适应这一变化，引入了状态变量法，也称内部法。
因此，系统的性能必须在设计、制造和使用环节应用先进的可靠性工程技术来提高，在多状态、多阶段理论框架下开展对其可靠性理论与方法的研究就显得尤为重要[2]。
本课题正是选取高阶典型的连续系统，首先用基本元件仿真，传输函数，状态空间的三种方法对系统建模，设置相应参数，改变输入信号类型，幅度和相位，采样间隔等，观察输出波形，并且编写M文件程序，调用模型分析系统频谱，使信号的各种性能能很好地得到研究。
1状态空间分析方法的概述
状态空间法是一种基于解答空间的问题表示和求解方法，它是以状态和操作符为基础的。
在利用状态空间图表示时，从某个初始状态开始，每次增加一个操作符，递增地建立起操作符的试验序列，直到达到目标状态为止。
由于状态空间法需要扩展过多的节点，容易出现“组合爆炸”，因而只适用于表示比较简单的问题。
系统的实现:
根据系统的外部描述构造一个内部结构，要求既保持外部描述的输入输出关系，又要将系统的内部结构确定下来。
这是一个复杂的问题，但也是一个非常重要的问题。
一方面，描述系统输入输出关系的微分方程或传递函数可以用实验的方法得到，我们可以从输入输出关系描述建立状态空间描述，这是建立状态空间描述的一条途径[3-5]。
一般描述为:
（1.1）
设状态变量选为
（1.2）
则（1.3）
由公式1.2和1.3有
因此，系统的状态方程：
（1.4）
输出方程为y=x1，由微分方程表达为矩阵形式，如下：
（1.5）
（1.6）
2快速创建LTI连续系统状态空间模型的方法
2.1创建LTI连续系统传递函数的方法
对于一般的低阶电路系统，例如RLC三阶电路（π电路），要得到它的系统函数，可以先由把三阶电路系统映射到S域，可以计算出其等效阻抗，再在复频域进行求解，这样比较方便简单。
具体做法如下：
将电路映射到S域，在频域的电路如1图所示：
图1三阶RLC复频域电路
利用拉普拉斯变换进行分析：
电阻不变，电感L变为Ls，电容C变为1/Cs,由复频域电路图建立复频域代数方程:
（2.1）
其等效阻抗：
（2.2）
输出像函数：
（2.3）
代入R,L,C的值，最后的输出象函数：
（2.4）
由输出象函数得系统函数：
（2.5）
2.2构造描述LTI连续系统的状态空间模型矩阵
要由传递函数求解状态空间矩阵的方法，有两种。
一种是通过迭代法演算，先将传递函数写成微分（差分）方程，再根据定义或画出流程图进行计算，这种方法只适合低阶系统，高阶系统会变的复杂，繁琐，所以应用不广泛，但不抽象，有利于理解。
另一种是通过Matlab编程直接求解，它简单，快捷，但具体算法不明确[6]。
若已求得H（S）的分子、分母系数分别为nu、de,在Matlab里可以调用函数tf2ss进行转换，tf2ss（nu,de）是一个函数，参数为nu、de，函数的返回值有4个[7]；如二阶连续系统的系统函数为H（s）=（d*s+c）/（s^2+a*s+b）;nu=[dc];de=[1ab];转换为状态空间模型；x=A*x+B*u;y=C*x+D*u;其中的四个矩阵可由[A,B,C,D]=tf2ss（nu,de）；给出对应初值条件；具体做法如下：
设de=[1.0000,-4.3244,5.5678,-3.3757,2.7467,-2.3446,1.8653,-1.4578];
nu=[-0.0000,0.0003,0.0053,0.0332,0.0384,0.0193,0.0000,0];
运行结果如下：
A=
4.3244-5.56783.3757-2.74672.3446-1.86531.4578
1.0000000000
01.000000000
001.00000000
0001.0000000
00001.000000
000001.00000
B=
1
0
0
0
0
0
0
C=
0.00030.00530.03320.03840.019300
D=0
由此，便可由传递函数得到状态空间矩阵。
3用Simulink状态空间建模求解LTI系统数值解的思路
3.1LTI连续系统的描述
时间和各个组成部分的变量都具有连续变化形式的系统称为连续系统。
连续系统的运动过程常可用一个或一组反映其变量间的相互联系和因果关系的微分方程来描述，通常称这一方程为它的数学模型。
当微分方程的系数为常数时称为定常系统，当系数随时间而变化时则称为时变系统。
如果连续系统的数学模型具有常微分方程的形式，则该系统属于集中参数系统；而当具有偏微分方程的形式时，相应的连续系统属于分布参数系统[8]。
3.2创建系统的Simulink状态空间模型
在仿真过程中，开始也学习认识到Matlab/Simulink系统建模设计与仿真基础知识，例如Matlab在系统建模中的应用有求取系统传递函数（矩阵），求取幅频特性和相频特性，系统输出响应，求系统状态空间模型，数值分析与函数逼近等等[9]。
数学仿真采用数学模型，用数学语言对系统的特性进行描述，其工作过程是：
（1）建立系统的数学模型；
（2）建立系统仿真模型，即设计算法，并转化为计算机程序，使系统的数学模型能为计算机所接受并能在计算机上运行；
（3）运行仿真模型，进行仿真试验，再根据仿真试验的结果进一步修正系统的数学模型和仿真模型。
模型基本结构：
一个典型的Simulink模型包括如下三种类型的元素：
信号源模块，被模拟的系统模块，输出显示模块。
两种Simulink运行仿真的方法使用窗口运行仿真；使用Matlab命令运行仿真。
使用窗口运行仿真优点：
人机交互性强，不必记住繁琐的命令语句即可进行操作。
使用窗口运行仿真主要可以完成以下一些操作。
（1）设置仿真参数
（2）应用仿真参数
（3）启动仿真
（4）停止仿真
（5）中断仿真
（6）仿真诊断
状态空间模型是动态时域模型，以隐含着的时间为自变量。
状态空间模型在经济时间序列分析中的应用正在迅速增加[10]。
根据以上准备工作，先在Simulink中建模如下，信源选择方波。
如图2所示：
图2用状态空间创建仿真模块
3.3模块内部参数设置及数据存储
在模块创建中，各个模块的内部参数设置如下。
信号源模块是方波信号，内部设置参数如图3。
状态空间内部参数设置为矩阵A,B,C,D，初始条件为x0.如图4。
数据存储即工作空间如图5和图6所示，其中时间为t为结算时间，x表示states，是指输入变量的存储，y表示output,是指输出变量的存储。
图3方信号源内部参数设置
图4状态空间模型内部参数设置
图5数据存储
图6数据存储
4利用Simulink状态空间建模求解LTI系统的优缺点
4.1状态空间建模求解LTI系统的优点
在科学研究和工程应用中更多地需要系统的数值解。
调用Matlab的Simulink工具包，可以用模块图标方法来编程，并通过运行系统模型文件的方法直接给出连续系统的数值解，而不需要将连续系统转换为离散系统再求解，省时省力。
其优点为：
可应用范围广，实用性强，可刻画系统的内部特征，用“一阶微分方程组”来描述系统的动态特性。
它的主要特点是：
（1）利用描述系统内部特性的状态变量替代了仅能描述系统外部特性的系统函数，能完整地揭示系统的内部特性，从而使得控制系统的分析和设计产生根本性的变革。
（2）便于处理多输入—多输出系统。
（3）一阶微分（或差分）方程组便于计算机数值计算。
4.2状态空间建模求解LTI系统的缺点
（1）不直观，几何、物理意义不明显：
只能给出系统的数值解，没有办法直接给出系统的解析解，只能根据系统的数值解及频谱分析给出近似的解析解。
（2）对数学模型要求很高：
而实际中往往难以获得高精度的模型，这妨碍了它推广和应用。
（3）一阶微分（或差分）方程组便于计算机数值计算。
（4）容易推广用于时变系统和非线性系统。
5连续系统Simulink状态空间建模分析方法程序设计的思路
5.1调用模型文件及编程求解系统响应
求解响应Matlab程序如下：
%二阶连续系统用Simulink状态空间建模，可分别求0输入和0状态响应
%y''（t）+a*y'（t）+b*y（t）=c*f（t）+d*f'（t）,y（0）=0,y'（0）=0
%A=[-a-b;10];B=[1;0];C=[dc];D=[0];x0=[q;p];
%x'（t）=A.x（t）+Bf[t];x（0）=x0;y（t）=C.x（t）;
%6个常数a,b,c,d,p,q完全描述了二阶连续系统
clc;clearall;closeall;
%方波信源的幅度u,周期T，高电平宽度百分比数hw
u=5;T=10;hw=5;
%二阶连续系统结构数据a,b,c,d及初值数据p,q
a=1.6;b=8*pi;c=1;d=4;p=2;q=5;
%定步长求解，步长为Ts,求解时区为Tsp
Ts=0.01;tf=20;L=round（tf/Ts）;Tsp=[0,（L-1）*Ts];
%=============0输入响应yzi（t）================
u=0;A=[-a-b;10];B=[0;0];C=[01];D=[0];x0=[q;p];
[t,x,y]=sim（'jxl00',Tsp）;
xin=y（:
1）;%输入信号采样放在y的第一列
yzi=y（:
2）;%0输入响应yzi（t）放在y的第二列
%==============0状态响应yzs（t）================
u=5;A=[-a-b;10];B=[1;0];C=[dc];D=[0];x0=[0;0];
[t,x,y]=sim（'jxl00',Tsp）;
yzs=y（:
2）;%0状态响应yzi（t）放在y的第二列
%==============全响应y（t）====================
yt=yzi+yzs;
%===============绘制波形======================
subplot（311）;plot（t,yzi,'b-'）;xlabel（'t'）;ylabel（'yzi（t）'）;
title（'二阶连续系统在周期方波激励下的0输入响应'）;
subplot（312）;plot（t,yzs,'r-'）;xlabel（'t'）;ylabel（'yzs（t）'）;
title（'二阶连续系统在周期方波激励下的0状态响应'）;
subplot（313）;plot（t,yt,'g-'）;xlabel（'t'）;ylabel（'y（t）'）;
title（'二阶连续系统在周期方波激励下的全响应'）;
程序运行结果如图7所示：
图7高阶系统在激励下的各种响应
5.2分析系统的频谱与相位
程序如下：
%===============频谱分析======================
fs=1/Ts;tao=Tspan
（2）+Ts;df=1/tao;
L=length（t）;L2=round（L/2）;
Y1=fft（yzi）*（2/L）;Y1
（1）=Y1
（1）/2;AY1=abs（Y1（1:
L2））;Am1=max（AY1）;
BY1=angle（Y1（1:
L2））;Bm1=max（BY1）;
Y2=fft（yzs）*（2/L）;Y2
（1）=Y2
（1）/2;AY2=abs（Y2（1:
L2））;Am2=max（AY2）;
BY2=angle（Y2（1:
L2））;Bm2=max（BY2）;
F=（0:
（L2-1））*df;
figure
（2）;
subplot（211）;stem（F,AY1,'r.'）;xlabel（'f/Hz'）;ylabel（'A1（f）'）;title（'系统零输入响应频谱'）;
subplot（212）;stem（F,AY2,'b.'）;xlabel（'f/Hz'）;ylabel（'A2（f）'）;title（'系统零状态响应频谱'）;
figure（3）;
subplot（211）;stem（F,BY1,'r.'）;xlabel（'f/Hz'）;ylabel（'B1（f）'）;title（'系统零输入响应相位'）;
subplot（212）;stem（F,BY2,'g.'）;xlabel（'f/Hz'）;ylabel（'B1（f）'）;title（'系统零状态响应相位'）;
频谱分析如图8所示：
图8高阶系统的频谱分析
6状态空间分析方法的应用实例
6.1实际连续系统的描述
设我们要分析如图9所示的二阶RLC电路，
图9二阶RLC电路
已知：
L=1H,C=2F,R=1Ω，求该系统函数。
解：
其中各元件的伏安关系为
，
（6.1）
以
为响应,则得到以
为响应的微分方程为：
（6.2）
将电阻电路映射到S域，利用拉普拉斯变换进行分析：
电阻不变，电感电感L变为
，电容C变为
，则求出系统函数：
（6.3）
代入数据，即可得
6.2在程序中设置参数
系统函数的分子nu=[000.5];分母de=[10.50.5];调用函数tf2ss进行转换，[A,B,C,D]=tf2ss（nu,de）；即可求出该函数的状态空间矩阵。
部分程序如下：
clc;clearall;closeall;
a=0.5;b=0.5;c=0;d=0.5;
nu=[cde];
de=[1ab];
[A,B,C,D]=tf2ss（nu,de）
运行结果为：
A=[-0.5-0.5;10]
B=[1;0]
C=[00.5]
D=[0]
6.3运行程序，求解系统
编写完整程序，运行程序，求出系统的数值解，并且描述该数值解。
clc;clearall;closeall;
%方波波波信源的幅度u,周期T
u=5;T=10;hw=5;
%二阶连续系统结构数据a,b,c,d及初值数据p,q
a=0.5;b=0.5;c=0;d=0.5;p=2;q=5;
%定步长求解，步长为Ts,求解时区为Tsp
Ts=0.01;tf=20;L=round（tf/Ts）;Tsp=[0,（L-1）*Ts];
%=============0输入响应yzi（t）================
u=0;A=[-a-b;10];B=[0;0];C=[01];D=[0];x0=[q;p];
[t,x,y]=sim（'jxl00',Tsp）;
xin=y（:
1）;%输入信号采样放在y的第一列
yzi=y（:
2）;%0输入响应yzi（t）放在y的第二列
%==============0状态响应yzs（t）================
u=5;A=[-a-b;10];B=[1;0];C=[dc];D=[0];x0=[0;0];
[t,x,y]=sim（'jxl00',Tsp）;
yzs=y（:
2）;%0状态响应yzi（t）放在y的第二列
%==============全响应y（t）====================
yt=yzi+yzs;
%===============画波形======================
subplot（311）;plot（t,yzi,'b-'）;xlabel（'t'）;ylabel（'yzi（t）'）;
title（'二阶连续系统在周期方波激励下的0输入响应'）;
subplot（312）;plot（t,yzs,'r-'）;xlabel（'t'）;ylabel（'yzs（t）'）;
title（'二阶连续系统在周期方波激励下的0状态响应'）;
subplot（313）;plot（t,yt,'g-'）;xlabel（'t'）;ylabel（'y（t）'）;
title（'二阶连续系统在周期方波激励下的全响应'）;
运行结果如图10和图11所示：
图10二阶系统在激励下的各种响应图11高阶系统的频谱分析
列举部分y（t）在t时间的数值如表1所示：
表1系统在部分t时刻的函数值
时间t
0.01
0.02
0.03
0.04
0.05
0.06
0.07
0.08
0.09
0.01
……
y（t）
5.3
5.1
3.2
0.6
-1.0
-1.7
-1.5
-1.0
-0.2
0.2
……
结语
现代社会是信息的社会，与人们生活密切相关的都是通信领域，比如决定现代社会发展方向的电脑以及几乎每个同学都拥有的手机等，充分说明了通信系统的重要性。
学会Matlab软件，并对通信系统应用和开发的设计思想有了更进一步的了解和掌握。
通过这次毕业设计我发现了毕业设计不仅仅是对以前所学知识的检验，也是对自己能力的提高。
这次毕业设计使我明白自己很多方面的知识还是比较欠缺，自己要学习的东西还有很多，以前老是有点眼高手低。
通过这次毕业设计，我明白了学习是长期积累的过程，无论是在以后的工作还是生活中都应该不断的学习，努力提高自己的知识和综合素质。
在设计的过程中我总会遇到这样那样的问题，我有失落过，烦恼过，悲伤过，但我明白这又是我人生中的一大挑战，角色的转换，这除了有较强的适应力和乐观的生活态度外，更重要
升级会员 