开放实验系统建模.docx
《开放实验系统建模.docx》由会员分享，可在线阅读，更多相关《开放实验系统建模.docx（33页珍藏版）》请在冰点文库上搜索。
开放实验系统建模
实验一 LTI系统建模与仿真
微分方程的符号数学求解
一、 实验目的
1. 熟悉连续时间系统的单位冲激响应、阶跃响应的意义及求解方法
2. 熟悉连续（离散）时间系统在任意信号激励下响应的求解方法
3. 熟悉应用MATLAB实现求解微分方程的方法
二、 实验原理
1.连续时间系统
对于连续的LTI系统，当系统输入为f（t），输出为y（t），则输入与输出之间满足如下的线性常系数微分方程：
，当系统输入为单位冲激信号δ（t）时产生的零状态响应称为系统的单位冲激响应，用h（t）表示。
若输入为单位阶跃信号ε（t）时，系统产生的零状态响应则称为系统的单位阶跃响应，记为g（t）。
二阶电路微分方程求解（符号数学求解）
uc（0）=3,i（0）=0
解：
零输入响应：
当L=0.5;C=0.02，本二阶电路的谐振频率f=1.5915，周期T=0.6283
>>uc=dsolve（'L*C*D2u+R*C*Du+u=Us','u（0）=3','Du（0）=0'）
uc=
-1/2*exp（-1/2*（C*R-（C^2*R^2-4*C*L）^（1/2））/C/L*t）*（（C^2*R^2-4*C*L）^（1/2）*Us-3*（C^2*R^2-4*C*L）^（1/2）+C*R*Us-3*C*R）/（C^2*R^2-4*C*L）^（1/2）+1/2*exp（-1/2*（C*R+（C^2*R^2-4*C*L）^（1/2））/C/L*t）*（-（C^2*R^2-4*C*L）^（1/2）*Us+3*（C^2*R^2-4*C*L）^（1/2）+C*R*Us-3*C*R）/（C^2*R^2-4*C*L）^（1/2）+Us
>>L=0.5;C=0.02;R=0.3;Us=0;%R=10为临界状态
>>y=subs（uc）
y=
1/19982*exp（（-3/10+1/5000*（-9991）^（1/2）*250000^（1/2））*t）*（-3/250000*（-9991）^（1/2）*250000^（1/2）-9/500）*（-9991）^（1/2）*250000^（1/2）-1/19982*exp（（-3/10-1/5000*（-9991）^（1/2）*250000^（1/2））*t）*（3/250000*（-9991）^（1/2）*250000^（1/2）-9/500）*（-9991）^（1/2）*250000^（1/2）
>>ezplot（y,[0,5]）
冲击响应：
>>uc=dsolve（'L*C*D2u+R*C*Du+u=dirac（t）','u（0）=3','Du（0）=0'）
uc=
1/2*exp（-1/2*（R*C-（R^2*C^2-4*L*C）^（1/2））/L/C*t）*（-2*C+3*R*f*C+3*（R^2*C^2-4*L*C）^（1/2）*f）/（R^2*C^2-4*L*C）^（1/2）/f-1/2*exp（-1/2*（R*C+（R^2*C^2-4*L*C）^（1/2））/L/C*t）*（-2*C+3*R*f*C-3*（R^2*C^2-4*L*C）^（1/2）*f）/（R^2*C^2-4*L*C）^（1/2）/f-C*heaviside（t）*（-exp（1/2*（R*C+（-C*（-R^2*C+4*L））^（1/2））/L/C*t）+exp（1/2*（R*C-（-C*（-R^2*C+4*L））^（1/2））/L/C*t））*exp（-t*R/L）/（-C*（-R^2*C+4*L））^（1/2）/f
>>L=0.5;C=0.02;R=0.3;f=1000;
>>y=subs（uc）
y=
-1/19982*exp（（-3/10+1/5000*（-9991）^（1/2）*250000^（1/2））*t）*（-1/25+9/500/t+3/250000*（-9991）^（1/2）*250000^（1/2）/t）*（-9991）^（1/2）*250000^（1/2）*t+1/19982*exp（（-3/10-1/5000*（-9991）^（1/2）*250000^（1/2））*t）*（-1/25+9/500/t-3/250000*（-9991）^（1/2）*250000^（1/2）/t）*（-9991）^（1/2）*250000^（1/2）*t+1/499550*heaviside（t）*（-exp（（3/10+1/5000*（-9991）^（1/2）*250000^（1/2））*t）+exp（（3/10-1/5000*（-9991）^（1/2）*250000^（1/2））*t））*exp（-3/5*t）*（-9991）^（1/2）*250000^（1/2）*t
>>ezplot（y）
阶跃响应：
>>uc=dsolve（'L*C*D2u+R*C*Du+u=heaviside（t）','u（0）=3','Du（0）=0'）
uc=
3/2*exp（-1/2*（R*C-（R^2*C^2-4*L*C）^（1/2））/L/C*t）*（R*C+（R^2*C^2-4*L*C）^（1/2））/（R^2*C^2-4*L*C）^（1/2）-3/2*exp（-1/2*（R*C+（R^2*C^2-4*L*C）^（1/2））/L/C*t）*（R*C-（R^2*C^2-4*L*C）^（1/2））/（R^2*C^2-4*L*C）^（1/2）-1/2*（（-R*C+（-C*（-R^2*C+4*L））^（1/2））*exp（1/2*（R*C-（-C*（-R^2*C+4*L））^（1/2））/L/C*t）+（R*C+（-C*（-R^2*C+4*L））^（1/2））*exp（1/2*（R*C+（-C*（-R^2*C+4*L））^（1/2））/L/C*t）-2*（-C*（-R^2*C+4*L））^（1/2）*exp（t*R/L））/（-C*（-R^2*C+4*L））^（1/2）*heaviside（t）*C*exp（-t*R/L）/f
信号源Us=3*t：
uc=dsolve（'L*C*D2u+R*C*Du+u=3*t','u（0）=3','Du（0）=0'）
L=0.5;C=0.02;R=0.03;f=1000;
y=subs（uc）
ezplot（y,[05]）
gridon
信号源Us=exp（j*3*t）
uc=dsolve（'L*C*D2u+R*C*Du+u=exp（j*3*t）','u（0）=3','Du（0）=0'）
L=0.5;C=0.02;R=0.03;
y=subs（uc）
ezplot（abs（y）,[05]）
gridon
三、 实验内容
1.已知描述系统的微分方程和激励信号e（t）分别如下，试用解析方法求系统的单位冲激响应h（t）和零状态响应r（t），并用MATLAB绘出系统单位冲激响应和系统零状态响应的波形，验证结果是否相同。
①
；
②
；
③
；
【答
①y=dsolve（'D2y+4*Dy+4*y=-exp（-t）+dirac（t）','Dy（0）=0','y（0）=0'） ;ezplot（y） ;
②y=dsolve（'D2y+2*Dy+26*y=dirac（t）','Dy（0）=0','y（0）=0'）;ezplot（y） ;
③y=dsolve（'D2y+4*Dy+3*y=exp（-2*t）*heaviside（t）','Dy（0）=0','y（0）=0'）;ezplot（y） ;
】
2.如下图所示的电路中，已知
，
，且两电感上初始电流分别为
，如果以电阻
上电压
作为系统输出，请求出系统在激励
（v）作用下的全响应。
四、 预习要求
1.熟悉系统响应的求解方法。
2.了解MATLAB语言中关于符号数学运算及绘图。
五、实验报告要求
1.建立系统模型，理论计算微分方程的齐次解，特解，全响应的表达式。
由初始条件确定出积分常量。
并写出解题过程。
2.记录仿真结果（包括数据和波形）。
3.写出程序清单。
4.实验总结（收获及体会）
实验二LTI系统建模及时域仿真（特征根法求解）
-------零输入响应与零状态响应
掌握信号经过LTI系统的时域分析方法。
巩固已经学过的知识，加深对知识的理解和应用，加强学科间的横向联系，学会应用MATLAB对实际问题进行仿真。
学会对带有非零起始状态的LTI系统进行仿真。
二、设计要求
（1）根据实际问题建立系统的数学模型。
电路如下，建立系统的数学模型，并计算其完全响应；
（2）用MATLAB描述此系统；
（3）仿真实现并绘制输出信号的波形。
仿真结果与理论值比较。
三、设计方法与步骤：
一般的连续时间系统分析有以下几个步骤:
①求解系统的零输入响应;②求解系统的零状态响应;③求解系统的全响应;④分析系统的卷积；⑤画出它们的图形.下面以具体的微分方程为例说明利用MATLAB软件分析系统的具体方法.
1．连续时间系统的零输入响应的数值计算（特征根法）
描述n阶线性时不变（LTI）连续系统的微分方程为：
已知y及各阶导数的初始值为y（0）,y
（1）（0），…，y（n-1）（0），求系统的零输入响应。
建模：
当LIT系统的输入为零时，其零输入响应为微分方程的齐次解（即令微分方程的等号右端为零），其形式为（设特征根均为单根）
其中p1,p2,…,pn是特征方程a1λn+a2λn-1+…+anλ+an=0的根，它们可以用root（a）语句求得。
各系数Ci由y及其各阶导数的初始值来确定。
对此有
………………………
写成矩阵形式为:
P1n-1C1+P2n-1C2+…+Pnn-1Cn=Dn-1y0
即V•C=Y0
其解为：
C=V\Y0
式中
V为范德蒙矩阵，在matlab的特殊矩阵库中有vander。
以下面式子为例：
初始条件为
；
MATLAB程序：
a=input（'输入分母系数（微分方程左边系数）：
a=[a1,a2,...]='）;
n=length（a）-1;
Y0=input（'输入初始条件向量Y0=[y0,Dy0,D2y0,...]='）;
p=roots（a）;%求特征根
V=rot90（vander（p））;
c=V\Y0';%c为积分常量
dt=input（'微分间隔dt='）;
te=input（'时间终止值te='）;
t=0:
dt:
te;
y=zeros（1,length（t））;
fork=1:
n
y=y+c（k）*exp（p（k）*t）;
end
plot（t,real（y））;grid
xlabel（'t'）;ylabel（'y'）;
title（'零输入响应'）;
例如：
用这个通用程序来解一个二阶系统，运行此程序并输入
a=[1,2,3]Y0=[2,7]dt=0.002te=7
结果如下图所示。
根据图可以分析零输入响应，它的起始值与输入函数无关，只与它的初始状态值有关，其起始值等于y（0_）的值。
随着时间的推移，最后零输入响应的值无限的趋近于0。
2．连续时间系统零状态响应的数值计算（求特征根法）
我们知道，LTI连续系统可用如下所示的线性常系数微分方程来描述，
例如，对于以下方程:
可用
输入函数
，得出它的冲激响应h，再根据LTI系统的零状态响应y（t）是激励u（t）与冲激响应h（t）的卷积积分。
注意，如果微分方程的左端或右端表达式中有缺项，则其向量a或b中的对应元素应为零，不能省略不写，否则出错。
求系统的零状态响应
及初始状态
。
输入函数
。
建模
先求出系统的冲激响应，写出其特征方程
求出其特征根为p1和p2，及相应的留数r1，r2；则冲激响应为
输出y（t）可用输入u（t）与冲激响应h（t）的卷积求得。
MATLAB程序：
a=input（'输入分母系数a=[a1,a2,...]='）;
b=input（'输入输入信号系数b=[b1,b2,...]='）;
dt=input（'采样时间间隔：
dt='）;
te=input（'终止时间：
te='）;%te为终止时间
t=0:
dt:
te;
u=input（'输入信号函数：
u='）;
[r,p,k]=residue（b,a）;
h=r
（1）*exp（p
（1）*t）+r
（2）*exp（p
（1）*t）;%由传输函数求冲激响应,laplace逆变换
subplot（2,1,1）,plot（t,h）;grid
title（'冲激响应'）;
y=conv（u,h）*dt;
subplot（2,1,2）,
plot（t,y（1:
length（t）））;grid
title（'零状态响应'）;
程序运行结果
执行这个程序，取a=[1,2,3]b=[4,6]dt=0.001te=7
得出图形如下：
由于初始状态为零，所以零状态的起始值也为零，即h（t）包含了连续系统的固有特性，与系统的输入无关。
只要知道了系统的冲激响应，即可求得系统在不同输入时产生的输出。
因此，求解系统的冲激响应h对进行连续时间系统的分析具有非常重要的意义
3．连续时间系统的全响应计算
上面通过对LTI系统函数的描述，我们可以得知：
如果在系统的初始状态不为零，在激励f（t）的作用下，LTI系统的响应称为全响应，它是零输入响应和零状态响应之和，即
故可先求出零输入响应和零状态响应，再把两者相加，得到全响应。
但简单的相加可能由于零输入与零状态的矩阵不同而不能的出正确的结果，这就需要对矩阵进行截取，使它们的矩阵尺寸相同。
例如，对于以下方程：
初始值为：
y（0_）=2,
；
输入函数为：
求它的全响应。
建模
先根据零输入响应的求法，得出零输入响应y1（t）。
再根据零状态响应的求法，得出零状态响应y2（t）。
最后，全响应y等于零输入响应y1（t）加上零状态响应y2（t），得出全响应。
MATLAB程序：
a=input（'输入分母系数（微分方程左边系数）：
a=[a1,a2,...]='）;
n=length（a）-1;%n阶微分方程
Y0=input（'输入初始条件向量Y0=[y0,Dy0,D2y0,...]='）;
b=input（'输入分子系数（微分方程右边系数）：
b=[b1,b2,...]='）;
u=input（'输入信号函数u='）;
dt=input（'dt='）;te=input（'te='）;
t=0:
dt:
te;
p=roots（a）;V=rot90（vander（p））;c=V\Y0';
y1=zeros（1,length（t））;
fork=1:
ny1=y1+c（k）*exp（p（k）*t）;end
te=t（end）;
dt=te/（length（t）-1）;
[r,p,k]=residue（b,a）;
h=r
（1）*exp（p
（1）*t）+r
（2）*exp（p
（1）*t）;
y2=conv（u,h）*dt;
y=y1（1:
length（t））+y2（1:
length（t））;
figure
（1）;
subplot（3,1,1）,plot（t,y1）,grid
xlabel（'t'）;ylabel（'y1'）;title（'零输入响应'）;
subplot（3,1,2）,plot（t,y2（1:
length（t）））;grid
xlabel（'t'）;ylabel（'y2'）;title（'零状态响应'）;
subplot（3,1,3）,plot（t,y）,grid
xlabel（'t'）;ylabel（'y'）;title（'全响应响应'）;
程序运行结果
执行程序，取a=[1,2,3]Y0=[2,7]b=[1,4,6]
u=sin（2*t）+cos（5*t）dt=0.001te=7
结果如下图：
在零输入响应中任一时刻取值y1，在零状态响应的对应时刻取值y2，再在全响应的对应时刻取值y。
可以得出：
y=y1+y2。
三、 实验内容
1.已知描述系统的微分方程和激励信号e（t）分别如下，试用解析方法求系统的单位冲激响应h（t）和零状态响应r（t），并用MATLAB绘出系统单位冲激响应和系统零状态响应的波形，验证结果是否相同。
①
；
②
；
③
；
2.如下图所示的电路中，已知
，
，且两电感上初始电流分别为
，如果以电阻
上电压
作为系统输出，请求出系统在激励
（v）作用下的全响应。
四、 预习要求
1.熟悉系统响应的求解方法
2.了解MATLAB语言中关于系统分析的方法：
五、实验报告要求
1.理论上计算出系统的零输入响应、零状态响应、全响应的表达式，并写出解题过程。
2.记录仿真结果（包括数据和波形）。
3.写出程序清单。
4.实验总结（收获及体会）
实验三 LTI系统的系统函数及响应-单位冲激响应
一、 实验目的
1. 熟悉连续时间系统的单位冲激响应、阶跃响应的意义及求解方法
2. 熟悉连续（离散）时间系统在任意信号激励下响应的求解方法
3. 熟悉应用MATLAB实现求解系统响应的方法
二、 实验原理
1.连续时间系统
对于连续的LTI系统，当系统输入为f（t），输出为y（t），则输入与输出之间满足如下的线性常系数微分方程：
，当系统输入为单位冲激信号δ（t）时产生的零状态响应称为系统的单位冲激响应，用h（t）表示。
若输入为单位阶跃信号ε（t）时，系统产生的零状态响应则称为系统的单位阶跃响应，记为g（t），如下图所示。
系统的单位冲激响应h（t）包含了系统的固有特性，它是由系统本身的结构及参数所决定的，与系统的输入无关。
我们只要知道了系统的冲激响应，即可求得系统在不同激励下产生的响应。
因此，求解系统的冲激响应h（t）对我们进行连续系统的分析具有非常重要的意义。
在MATLAB中有专门用于求解连续系统冲激响应和阶跃响应，并绘制其时域波形的函数impulse（）和step（）。
如果系统输入为f（t），冲激响应为h（t），系统的零状态响应为y（t），则有：
。
若已知系统的输入信号及初始状态，我们便可以用微分方程的经典时域求解方法，求出系统的响应。
但是对于高阶系统，手工计算这一问题的过程非常困难和繁琐。
在MATLAB中，应用lsim（）函数很容易就能对上述微分方程所描述的系统的响应进行仿真，求出系统在任意激励信号作用下的响应。
lsim（）函数不仅能够求出连续系统在指定的任意时间范围内系统响应的数值解，而且还能同时绘制出系统响应的时域波形图。
以上各函数的调用格式如下：
⑴impulse（）函数：
（数值解。
相当于离散系统的impz函数）
函数impulse（）将绘制出由向量a和b所表示的连续系统在指定时间范围内的单位冲激响应h（t）的时域波形图，并能求出指定时间范围内冲激响应的数值解。
impulse（b,a） 以默认方式绘出由向量a和b所定义的连续系统的冲激响应的时域波形。
impulse（b,a,t0） 绘出由向量a和b所定义的连续系统在0~t0时间范围内冲激响应的时域波形。
impulse（b,a,t1:
p:
t2） 绘出由向量a和b所定义的连续系统在t1~t2时间范围内，并且以时间间隔p均匀取样的冲激响应的时域波形。
y=impulse（b,a,t1:
p:
t2） 只求出由向量a和b所定义的连续系统在t1~t2时间范围内，并且以时间间隔p均匀取样的冲激响应的数值解，但不绘出其相应波形。
⑵step（）函数（数值解）
函数step（）将绘制出由向量a和b所表示的连续系统的阶跃响应，在指定的时间范围内的波形图，并且求出数值解。
和impulse（）函数一样，step（）也有如下四种调用格式：
step（b,a）
step（b,a,t0）
step（b,a,t1:
p:
t2）
y=step（b,a,t1:
p:
t2）
上述调用格式的功能和impulse（）函数完全相同，所不同只是所绘制（求解）的是系统的阶跃响应g（t），而不是冲激响应h（t）。
⑶lsim（）函数：
（数值解。
相当于离散系统的filter函数）
根据系统有无初始状态，lsim（）函数有如下两种调用格式：
①系统无初态时，调用lsim（）函数可求出系统的零状态响应，其格式如下：
lsim（b,a,x,t） 绘出由向量a和b所定义的连续系统在输入为x和t所定义的信号时，系统零状态响应的时域仿真波形，且时间范围与输入信号相同。
其中x和t是表示输入信号的行向量，t为表示输入信号时间范围的向量，x则是输入信号对应于向量t所定义的时间点上的取样值。
y=lsim（b,a,x,t） 与前面的impulse和step函数类似，该调用格式并不绘制出系统的零状态响应曲线，而只是求出与向量t定义的时间范围相一致的系统零状态响应的数值解。
②系统有初始状态时，调用lsim（）函数可求出系统的全响应，格式如下：
lsim（A,B,C,D,e,t,X0） 绘出由系数矩阵A,B,C,D所定义的连续时间系统在输入为e和t所定义的信号时，系统输出函数的全响应的时域仿真波形。
t为表示输入信号时间范围的向量，e则是输入信号e（t）对应于向量t所定义的时间点上的取样值，X0表示系统状态变量X=[x1,x2,…..xn]'在t=0时刻的初值。
[Y,X]=lsim（A,B,C,D,e,t,X0） 不绘出全响应波形，而只是求出与向量t定义的时间范围相一致的系统输出向量Y的全响应以及状态变量X的数值解。
显然，函数lsim（）对系统响应进行仿真的效果取决于向量t的时间间隔的密集程度，t的取样时间间隔越小则响应曲线越光滑，仿真效果也越好。
说明：
（1）当系统有初始状态时，若使用lsim（）函数求系统的全响应，就要使用系统的状态空间描述法，即首先要根据系统给定的方式，写出描述系统的状态方程和输出方程。
假如系统原来给定的是微分方程或系统函数，则可用相变量法或对角线变量等方法写出系统的状态方程和输出方程。
（2）显然利用lsim（）函数不仅可以分析单输入单输出系统，还可以分析复杂的多输入多输出系统。
例题1：
若某连续系统的输入为e（t），输出为r（t），系统的微分方程为：
①求该系统的单位冲激响应h（t）及其单位阶跃响应g（t）。
②若
求出系统的零状态响应y（t）
分析：
①求冲激响应及阶跃响应的MATLAB程序：
a=[1 5 6];b=[3 2];
subplot（2,1,1）,impulse（b,a,4）
subplot（2,1,2）,step（b,a,4）
运行结果如右：
②求零状态响应的MATLAB程序：
a=[1 5 6];b=[3 2];
p1=0.01; %定义取样
升级会员 