第5章控制系统的数学建模与仿真

控制系统的数学模型关系到对系统各方面性能的仿真分析，建立控制系统的数学模型是分析和设计控制系统的首要工作。本章简要介绍时域建模方法、频域建模方法和神经网络建模方法，并通过仿真示例重点介绍仿真方法的使用。
数学模型有动态模型与静态模型之分。描述系统动态过程的方程式，如微分方程、偏微分方程、差分方程等，称为动态模型； 在静态条件下，即变量的各阶导数为零，描述系统各变量之间关系的方程式称为静态模型。
同一个物理系统，可以用不同的数学模型来表达。例如，实际的物理系统一般含有非线性特性，因此系统的数学模型就应该是非线性的。严格地讲，实际系统的参数不可能是集中的，所以系统的数学模型又应该用偏微分方程描述。但是非线性方程或偏微分方程对应模型的参数辨识相当困难，有时甚至不可能。因此，为了便于问题的求解，常常在误差允许的范围内，忽略次要因素，用简化的数学模型来表示实际的物理系统。这样，同一个系统就有完整、复杂的数学模型和简单、准确性较差的数学模型之分。一般情况下，在建立数学模型时，必须在模型的简化性与分析结果的精确性之间做出折中选择。
此外，数学模型的形式有多种。为了便于分析研究，可能某种形式的数学模型比另一种更合适。例如，在求解最优控制问题或多变量系统的问题时，采取状态变量表达式(即状态空间表达式)比较方便； 但是在对单输入、单输出系统的分析中，采用输入输出间的传递函数(或脉冲传递函数)作为系统的数学模型比较合适。
所以在建立系统数学模型时，必须： 
(1) 全面了解系统的特性，确定研究目的以及准确性要求，确定是否需要忽略一些次要因素而使系统数学模型简化，既不致造成数学处理上的困难，又不致影响分析的准确性。一般在条件允许时，最初尽可能采用简化的常系数线性数学模型，然后再在线性模型分析的基础上考虑被忽略因素所引起的误差，最后建立起系统比较完善准确的数学模型。但是必须指出，由于数学分析方法上的误差，复杂的数学模型不一定能带来预期的准确结果。
(2) 根据所应用的系统分析方法，建立相应形式的数学模型(微分方程、传递函数等)，有时还要考虑计算机仿真求解的要求。




建立系统的数学模型主要有两条途径。第一条途径是利用人们已有的关于系统的知识，采用机理的方法建立数学模型。机理法是一种推理方法，用这种方法建立模型时，是通过系统本身机理(物理、化学规律)的分析确定模型的结构和参数，从理论上推导出系统的数学模型。这种利用机理法得出的数学模型称为机理模型或解析模型。第二条途径是根据对系统的观察，通过测量所得到的大量输入、输出数据，推断出被研究系统的数学模型。这种方法称为实验法，利用实验法建立的数学模型称为经验模型。一般来讲，采用实验法建立的数学模型，是系统模型化问题的唯一解。而采用机理法时，能够满足观测到的输入、输出数据关系的系统模型却有很多个。





5．1时域建模方法及示例
本书介绍的时域建模方法，主要是采用机理法建模。连续系统的时域模型主要是指微分方程描述系统的模型，离散系统的时域模型主要是指差分方程描述系统的模型，以及相应的状态空间模型。
1． 微分方程建模
控制系统的运动状态可由微分方程式描述，而微分方程式就是系统的一种数学模型。关于建立系统微分方程式的一般步骤如下： 
(1) 在允许的条件下适当简化，忽略一些次要因素。
(2) 根据物理或化学定律，列出元件的原始方程式。这里所说的物理或化学定律，是指牛顿第二定律、能量守恒定律、物质不灭定律、克希霍夫定律等。
(3) 列出原始方程式中中间变量与其他因素的关系式。这种关系式可能是数学方程式，或是曲线图。它们在大多数场合是非线性的。若条件允许，则应进行线性化处理，否则按非线性系统对待。
(4) 将上述关系式代入原始方程式，消去中间变量，就得元件的输入、输出关系的方程式。若在步骤(3)不能进行线性化，则输入输出关系方程式将是比较复杂的非线性方程式。
(5) 同理，求出其他元件的方程式。
(6) 从所有元件的方程式中消去中间变量，最后的系统输入输出微分方程式。


图51弹簧质量阻尼器系统

【例51】如图51所示的一个弹簧质量阻尼器系统，当外力f(t)作用时，系统产生位移y(t)，要求写出系统在外力f(t)作用下的运动方程式。在此，f(t)是系统的输入，y(t)是系统的输出。
建模的过程如下： 
(1) 设运动部件质量用m表示，按集中参数处理，列出原始方程式。根据牛顿第二定律，有

f(t)－f1(t)－f2(t)=md2y(t)dt2(51)

式中，f1(t)——阻尼器阻力； 
f2(t)——弹簧弹力。
(2) f1(t)和f2(t)为中间变量，找出它们与其他因素的关系。由于阻尼器是一种产生黏性摩擦或阻尼的装置，活塞杆和缸体发生相对运动时，其阻力与运动方向相反，与运动速度成正比，故有

f1(t)=Bdy(t)dt(52)

式中，B——阻尼系数。
设弹簧为线性弹簧，则有

f2 (t)=Ky(t)(53)

式中，K——弹簧的弹性系数。
(3) 将式(52)和式(53)代入式(51)，经整理后就得描述该弹簧质量阻尼器系统的微分方程式为

md2y(t)dt2+Bdy(t)dt+ky(t)=f(t)(54)

式中，m、B、K均为常数，故上式为线性定常二阶微分方程，则此机械位移系统为线性定常系统。
一般列写微分方程式时，输出量及其各阶导数项列写在方程式左端，输入项列写在右端。由于一般物理系统均有质量、惯性或储能元件等物理可实现的原因，左端的导数阶次总比右端的高。在本例中，有质量M，又有吸收能量的阻尼器B，系统有两个时间常数，故左端导数项最高阶次为2。
【例52】建立图52所示的RLC串联电路输入、输出电压之间的微分方程。


图52RLC电路图

建模过程如下： 
(1) 确定系统的输入、输出量，输入端电压ui(t)为输入量，输出端电压uo(t)为输出量。
(2) 列写微分方程，设回路电流为i(t)，由基尔霍夫定律可得

uR+uL+uC=ui(55)

式中，uR(t)、uL(t)、uC(t)分别为R、L、C上的电压降。

又由

uR(t)＝Ri(t),uL(t)=Ldi(t)dt

可得

Ldi(t)dt+Ri(t)+uC(t)=ui(t)(56)


(3) 消去中间变量，得出系统的微分方程。考虑uC(t)=uo(t)，根据电容的特性可得

i(t)=CduC(t)dt=Cduo(t)dt(57)


将式(57)代入式(56)，可得到系统的微分方程为

LCd2uo(t)dt2+RCduo(t)dt+uo(t)=ui(t)(58)


比较式(54)和式(58)可知，当两个方程式的系数相同时，从动态性能角度来看，两个系统是相同的。这就有可能利用电气系统来模拟机械系统，进行试验研究。而且从系统理论来说，就有可能脱离系统的具体物理属性，进行普遍意义的研究。
一般而言，线性定常系统的微分方程所建模型可由下述n阶微分方程描述

andny(t)dtn+an－1dn－1y(t)dtn－1+…+a1dy(t)dt+a0y(t)

=bmdmu(t)dtm+bm－1dm－1u(t)dtm－1+…+b1du(t)dt+b0u(t)

式中，y(t)——系统输出量； 
u(t)——系统输入量； 
a0，a1，…，an和b0，b1，…，bm——与系统结构参数有关的常数。
2． 差分方程建模
由于数字计算机、微处理器的迅速发展和广泛应用，数字控制器在许多场合都取代了模拟控制器。由于数字控制器接收、处理和传送的都是数字信号，所以如果在控制系统中有一处或几处信号不是时间的连续函数，而是以离散的脉冲序列或数字脉冲序列形式出现，那么这样的系统则称为离散控制系统。通常将系统中的离散信号是脉冲序列形式的离散系统称为采样控制系统或脉冲控制系统； 将系统中的离散信号是数字序列形式的离散系统称为数字控制系统或计算机控制系统。
离散控制系统和连续控制系统相比，既有本质上的不同，又有分析研究方面的相似性。利用Z变换法研究离散系统，可以把连续系统中的一些概念和方法推广到线性离散系统的分析和设计中。
描述离散控制系统的时域数学模型主要是差分方程，这里通过例子简要介绍差分方程的基本概念及其建立方法。
与线性定常连续系统用线性微分方程描述相似，线性定常离散系统常用差分方程来描述。下面举例进行说明。
设离散控制系统的结构图如图53所示。


图53采样控制系统的结构图


在第k个采样时间间隔中，零阶保持器的输出为eh(t)=e(kT)，kT≤t≤(k+1)T。在该周期内的输出c(t)则为

c(t)=c(kT)+e(kT)(t－kT)(59)

式中，kT≤t≤(k+1)T。
根据c［(k+1)T］=c(kT)+Te(kT)，可简写为

c(k+1)=c(k)+Te(k)(510)


考虑到e(k)=r(k)－c(k)，将之代入式(510)，得

c(k+1)+(T－1)c(k)=Tr(k)(511)


上式就是如图53所示采样控制系统的差分方程模型。
更一般的描述，线性定常离散系统通常可用如下n阶前向差分方程来描述： 

c(k+n)+a1c(k+n－1)+…+an－1c(k+1)+anc(k)

=b0r(k+m)+b1r(k+m－1)+…+bm－1r(k+1)+bmr(k)

式中，n为系统阶次； k为第k个采样周期，a1，a2，…，an和b0，b1，…，bm是与系统结构参数有关的常数。
3． 状态空间表达建模
从20世纪50年代后期开始，Bellman等提出了状态变量法。在用状态空间法分析系统时，系统的动态特性是由状态变量构成的一阶微分方程组来描述的。在计算机上求解一阶微分方程组比求解与之相应的高阶微分方程容易得多，而且可以同时得到系统全部独立变量的响应，因而可以同时确定系统的全部内部运动状态。此外，状态空间法还可以方便地处理初始条件，可以用来分析设计多变量、时变和非线性控制系统，也可以应用于随机过程和采样数据系统，因此它是研究复杂控制系统的理论基础。现代控制理论是在引入状态和状态空间概念的基础上发展起来的。确定控制系统状态空间的描述，即建立状态空间的数学模型是一个基本问题。状态空间数学模型的建立需要理解如下基本概念。
状态： 动态系统的状态是指系统的过去、现在和将来的运动状态。状态需要一组必要而充分的数据来说明。如RLC电路所示系统，这个系统的状态就是RLC电路每一时刻的回路电流和输出电压。
状态变量： 系统的状态变量就是指可以完全确定系统运动状态的最小一组变量。一个用 n阶微分方程描述的系统，就有n个独立变量，求得这n个独立变量的时间响应，系统的运动状态也就被完全确立了。因此，系统的状态变量就是n阶系统的n个独立变量。
设x1(t)，x2(t)，…，xn(t)为系统的一组状态变量，满足下列两个条件： 在初始时刻t=t0，这组变量x1(t0)，x2(t0)，…，xn(t0)的值可以表示系统在初始时刻的状态； 当系统在t≥t0的输入和上述初始状态确定以后，状态变量便能完全确定系统在任何t≥t0时刻的行为。
对同一个系统来说，究竟选取哪些变量作为状态变量不是唯一的，关键是这些状态变量是相互独立的，且其个数应等于微分方程的阶数，因为微分方程的阶数唯一地取决于系统中独立储能元件的个数，因此，状态变量的个数就应等于系统独立储能元件的个数。还应该指出，状态变量不一定是物理上可量测或可观测的量，但通常选择易于测量或观测的量作为状态变量，这样在系统实现最佳控制规律时，就可以得到所有需要反馈的状态变量。
状态向量： 如果完全描述一个系统的动态行为需要n个状态变量x1(t)，x2(t)，…，xn(t)，那么这n个状态变量作分量所构成的向量X(t)叫作该系统的状态向量，记作X(t)=［x1(t)x2(t)…xn(t)］T。
状态空间： 以状态变量x1(t)，x2(t)，…，xn(t)为坐标所构成的n维空间称为状态空间。任何状态都可以用状态空间中的一个点来表示。即在特定时刻t，状态向量X(t)在状态空间中是一点。已知初始时刻t0的状态X(t0)，就得到状态空间中的一个初始点。随着时间的推移，状态X(t)将在状态空间中描绘出一条轨迹，称为状态轨迹。这一轨线的形状完全由系统在t0时刻的初始状态X(t0)和t≥t0的输入以及系统的动态结构唯一决定。状态向量的状态空间模型联系了向量的代数结构和几何结构。
状态方程： 描述系统状态变量与系统输入之间关系的一阶微分方程组称为状态方程。
对于式(58)所描述的RLC电路系统，若令x1(t)=u0(t)，x2(t)=x·1(t)，即取x1、x2为此系统的一组状态变量，则可以得到一阶微分方程组： 

x·1
x·2=01
－1LC－RLx1
x2+0
1LCui(512)


式(512)即为该系统的状态方程，上式可简写成

X·=AX+Bui (513)

式中： 

A=01
－1LC－RL，B=0
1LC，X=x1
x2


输出方程： 描述系统的状态变量与输出变量关系的一组代数方程称为输出方程。
对于式(58)所描述的RLC电路系统，指定uo(t)为输出，则有

uo(t)=x1(t)(514)

写成

uo=CTX(515)

式中,

CT=10，X=x1
x2

式(515)即为该系统的输出方程。
状态方程和输出方程一起构成一个系统动态的完整描述，称为系统的状态空间模型表达式。如式(513)和式(515)就是式(58)所描述的RLC电路系统的状态空间模型表达式。
一般地，对于一个复杂系统，它可以有r个输入，m个输出，此时状态方程为

x·1=a11x1+a12x2+…+a1nxn+b11u1+…+b1rur

x·2=a21x1+a22x2+…+a2nxn+b21u1+…+b2rur



x·n=an1x1+an2x2+…+annxn+bn1u1+…+bnrur(516)


输出方程不仅是状态变量的组合，而且在特殊情况下可以有输入向量的直接传递，因而有以下一般形式： 

y1=c11x1+c12x2+…+c1nxn+d11u1+…+d1rur

y2=c21x1+c22x2+…+c2nxn+d21u1+…+d2rur



ym=am1x1+am2x2+…+amnxn+dm1u1+…+dmrur(517)


多输入多输出系统状态空间表达式的向量矩阵形式表示为

x·=Ax+Bu


y=Cx+Du(518)

式中，x和A分别为n维状态向量和n×n状态矩阵； u表示r维输入(或控制)向量； y表示m维输出向量； B表示n×r输入(或控制)矩阵； C表示m×n输出矩阵； D表示输入量的m×r直接传递矩阵。
下面简要叙述定常连续状态空间模型的离散化方法。
已知定常连续系统状态方程

x·=Ax+Bu(519)


在x(t0)以及输入u(t)作用下的解为

x(t)=Φ(t－t0)x(t0)+∫tt0Φ(t－τ)Bu(τ)dτ(520)


令t0=kT，有x(t0)=x(kT)=x(k)，令t=(k+1)T，则有
x(k+1)T=x(k+1)在t∈k,k+1时，有u(k)=u(k+1)为常数，于是其解为

x(k+1)=Φ(T)x(k)+∫(k+1)TkTΦ(k+1)T－τBdτu(k)(521)


记G(t)=∫(k+1)TkTΦ(k+1)T－τBdτ，为了便于计算G(T)，引入变换(k+1)T－τ=τ′，则有G(t)=∫0T－Φ(τ′)Bdτ′=∫T0Φ(τ)Bdτ，因此离散化系统的状态方程为

x(k+1)=Φ(T)x(k)+G(T)u(k)(522)

式中,Φ(T)根据连续系统的状态转移矩阵Φ(t)导出，即Φ(T)=Φ(t)|t=T。
离散化系统输出方程为

y(k)=Cx(k)+Du(k)(523)

4． 时域模型的仿真实现
(1) 对于微分方程模型的求解，第1章中已详细介绍了采用函数dsolve进行数值求解，采用函数solve进行解析求解的使用方法。
(2) 对于差分方程模型的求解，可以运用MATLAB提供的函数filter和filtic。这两个函数的使用说明如下。
函数filtic为函数filter求解差分方程设定初始条件。
filtic函数的使用格式为Z=filtic(B,A,Y,X)。其中B，A分别为差分方程 a(1)*y(n)+a(2)*y(n-1)+a(3)*y(n-2)+…+a(nb+1)*y(n-nb)=b(1)*x(n)+b(2)*x(n-1)+…+b(nb+1)*x(n-nb)的输入和输出变量前的系数； X，Y分别为输入和输出的初值。
函数filter用于求解差分方程。
filter函数的使用格式为Y=filter (B,A,X)。其中B，A分别为差分方程 a(1)*y(n)+a(2)*y(n-1)+a(3)*y(n-2)+…+a(nb+1)*y(n-nb)= b(1)*x(n)+b(2)*x(n-1)+…+b(nb+1)*x(n-nb)的输入和输出变量前的系数； X为输入序列。
【例53】求解差分方程y(n)-0．4y(n-1)-0．45y(n-2)=0．45x(n)+0．4x(n-1)-x(n-2)，初值为y(-1)=0，y(-2)=1，x(-1)=1，x(-2)=2。
求解程序如下： 

>> B=［0.45 0.4 -1］;

>> A=［1 -0.4 -0.45］;

>> x0=［1 2］;%设定初值

>> y0=［0 1］;

>> N=50;

>> n=［1:N-1］';

>> x=0.2.^n;%生成输入序列

>> Zi=filtic(B,A,y0,x0);%生成初始条件

>> ［y,Zf］=filter(B,A,x,Zi);

>> plot(n,x,'R-',n,y,'b--');

>> xlabel('n');ylabel('(n)--y(n)');

运行结果如图54所示。


图54例53程序运行结果


(3) 对于状态空间模型的求解，可以先用ss命令来建立状态空间模型，然后求解。
对于连续系统，其格式为 sys=ss(A,B,C,D)，其中A，B，C，D为描述线性连续系统的矩阵。当sys1是一个用传递函数表示的线性定常系统时，可以用命令sys=ss(sys1)将其转换成为状态空间形式。也可以用命令sys=ss(sys1,'min')计算出系统sys的最小实现。
和连续系统状态空间表达式的输入方法相类似，如果要输入离散系统的状态空间表达式： 

x(k+1)=Gx(k)+Hu(k)


y(k)=Cx(k)+du(k)(524)


首先需要输入矩阵G、H、C、d，然后输入语句sys=ss(G,H,C,d,T)，即可将其输入到MATLAB的工作空间中，并且用变量名来表示这个离散系统，其中T为采样时间。如果Gyu表示一个以脉冲传递函数描述的离散系统，也可以用ss(Gyu )命令，将脉冲传递函数模型转换成状态空间表达式。
在MATLAB中，还提供了函数c2d，其功能就是将连续时间的系统模型转换成离散时间的系统模型。其调用格式为sysd=c2d(sysc,T,method)。其中，输入参量sysc为连续时间的系统模型； T为采样周期； method用来指定离散化采用的方法。可以选用的离散化采用的方法有： 'zoh'为采用零阶保持器； 'foh'为采用一阶保持器； 'tustin'为采用双线性逼近方法； 'prewarm'为采用改进的tustin方法； 'matched'为采用SISO系统的零极点匹配方法； 当method默认，即调用格式为sysd=c2d(sysc,T)时，默认的方法是采用零阶保持器。
【例54】已知系统状态方程为
x·=01
－2－3x+0
1u，x(0)=1
0，求解系统在阶跃信号作用下的输出。
方法1，数值求解，求解程序如下： 

>>A=［0 1;-2 -3］;

>>B=［0 1］';

>> C=［1 0;0 1］;

>> D=［0 0］';

>>x0=［1 0］';

>> t=0:0.01:10;

>>u(length(t))=1;

>> u(:)=1;

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

>>lsim(sys,u,t,x0)

运行结果如图55所示。


图55例54程序运行结果


方法2，解析求解，求解程序如下： 

>>syms s t x0 x tao p p0;

>>A=［0 1;-2 -3］;

>>B=［0 1］';

>>I=［1 0;0 1］;

>>E=s*I-A;

>>C=det(E);

>>D=collect(inv(E));

>>p0=ilaplace(D);

>>x0=［1 0］';

>>x1=p0*x0;

>>p=subs(p0,'t',(t-tao));

>>F=p*B*1;

>>x2=int(F,tao,0,t);

>>x=collect(x1+x2)

运行结果如下： 

x =

2*exp(-t) - exp(-2*t) + (exp(-2*t)*(exp(t) - 1)^2)/2

2*exp(-2*t) - 2*exp(-t) + exp(-2*t)*(exp(t) - 1)

【例55】线性连续系统的状态方程为

x·=Ax+Bu


y=Cx+Du

其中，

A=010

001

-6-11-6
，B=10

2-1

02
，C=1-10

21-1
，D=00

00


采用零阶保持器将其离散化，设采样周期为0．1秒。求离散化的状态方程模型，及其在零初始条件下，该离散化系统的阶跃信号响应。
求解程序如下： 

>>A=［0 1 0;0 0 1;-6 -11 -6］;

>>B=［1 0;2 -1;0 2］;

>>C=［1 -1 0;2 1 -1］;

>>D=zeros(2);

>>T=0.1;

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

>>sysd=c2d(sys,T)

>>step(sysd)

运行结果如下： 

sysd =

a =

x1x2x3

x10.99910.09840.004097

x2-0.024580.95410.07382

x3-0.4429-0.83660.5112

b =

u1u2

x10.1099-0.004672

x20.1959-0.0902

x3-0.11640.1936

c =

x1x2x3

y11-10

y221-1

d =

u1u2

y100

y200

Sample time: 0.1 seconds

Discrete-time state-space model.

该离散化系统的阶跃信号响应如图56所示。


图56例55程序运行结果







5．2频域建模方法及示例
控制系统的复频域数学模型主要有传递函数、结构图、信号流图和频率特性，下面先简要叙述这些模型的基本概念，再介绍相应的MATLAB函数实现。
1． 传递函数模型
对于描述控制系统的微分方程模型，在给定初始条件下的情况下，可以通过求解微分方程直接得到系统的输出响应，但如果方程阶次较高，则计算很烦琐，从而给系统的分析设计带来不便。经典控制理论的主要研究方法，都不是直接利用求解微分方程的方法，而是采用与微分方程有关的传递函数模型进行描述。传递函数是经典控制理论中最重要的数学模型。利用传递函数不必求解微分方程就可研究初始条件为零的系统在输入信号作用下的动态性能。利用传递函数还可研究系统参数变化或结构变化对动态过程的影响，因而极大地简化了系统分析的过程。另外，还可以把对系统性能的要求转化为对系统传递函数的要求，使综合设计问题易于实现。


图57传递函数的图示

所谓传递函数，即线性定常系统在零初始条件下，系统输出量的拉普拉斯变换与输入量的拉普拉斯变换之比。如图57表示一个具有传递函数G(s)的线性系统，图中表明，系统输入量与输出量的关系可以用传递函数联系起来。
设线性定常系统为下述n阶线性常微分方程所描述： 

andnc(t)dtn+an－1dn－1c(t)dtn－1+…+a1dc(t)dt+a0

=bmdmr(t)dtm+bm－1dm－1r(t)dtm－1+…+b1dr(t)dt+b0(525)

式中，c(t)——系统输出量； r(t)——系统输入量。
在初始状态为零时，对上式两端取拉普拉斯变换，得

ansnC(s)+an－1sn－1C(s)+…+a1sC(s)+a0C(s)

=bmsmR(s)+bm－1sm－1R(s)+…+b1sR(s)+b0R(s)(526)


式(526)用传递函数可表述为

G(s)=C(s)R(s)=bmsm+bm－1sm－1+…+b1+b0ansn+an－1sn－1+…+a1+a0(527)

式中，C(s)——输出量的拉普拉斯变换； R(s)——输入量的拉普拉斯变换； G(s)——系统或环节的传递函数。
由于传递函数在经典控制理论中是非常重要的概念，故有必要对其性质、适用范围及表示形式等方面作出以下说明： 传递函数只适用于描述线性定常系统； 传递函数和微分方程一样，表征系统的运动特性，是系统的数学模型的一种表示形式，它和系统的运动方程是一一对应的。传递函数分子多项式系数及分母多项式系数，分别与相应微分方程的右端及左端微分算符多项式系数相对应。在零初始条件下，将微分方程的算符d/dt用复数s置换便得到传递函数； 反之，将传递函数多项式中的变量s用算法d/dt置换便得到微分方程； 传递函数是系统本身的一种属性，它只取决于系统的结构和参数，与输入量和输出量的大小和性质无关，也不反映系统内部的任何信息。且传递函数只反映系统的动态特性，而不反映系统物理性能上的差异，对于物理性质截然不同的系统，只要动态特性相同，它们的传递函数就具有相同的形式； 传递函数为复变量s的真有理分式，即n≥m，因为系统或元件总是具有惯性的，而且输入系统的能量也是有限的； 传递函数是在初始条件为零时定义的，因此，在非零初始条件下，传递函数不能完全表征系统的动态性能。另外，系统内部往往有多种变量，但传递函数只是通过系统输入量和输出量之间的关系来描述系统，而对内部其他变量的情况却无法得知。特别是某些变量不能由输出变量反映时，传递函数就不能正确表征系统的特征。现代控制理论采用状态空间法描述系统，引入了可控性和可观测性的概念，从而对控制系统进行全面的了解，可以弥补传递函数的不足； 传递函数G(s)的拉普拉斯逆变换是脉冲响应g(t)。
作为线性定常离散控制系统的数学模型，采用脉冲传递函数描述。脉冲传递函数的定义与连续系统的传递函数的定义类似。
以图58为例，如果系统的初始条件为零，输入信号为r(t)，采样后r*(t)的Z变换函数为R(z)，系统连续部分的输出为c(t)，采样后c*(t)的Z变换函数为C(z)，则线性定常离散系统的脉冲传递函数定义为系统输出采样信号的Z变换与输入采样信号的Z变换之比，记作

G(z)=C(z)R(z)(528)

此处，零初始条件是指t<0时，输入脉冲序列各采样值r(－T),r(－2T),…及输出脉冲序列各采样值c(－T),c(－2T),…均为零。


图58离散系统

由式（528）所描述的关系可以得知，输出的采样信号如下式所描述： 

c*(t)=Z－1［C(z)］=Z－1［G(z)R(z)］(529)


由于输入信号的Z变换R(z)通常为已知的，因此求c*(t)的关键在于求出系统的脉冲传递函数G(z)。
关于线性定常离散控制系统的脉冲传递函数，还需着重指出： G(s)是一个线性环节的传递函数，而G(z)表示的是线性环节与理想开关两者的组合的脉冲传递函数。如果不存在理想采样开关，那么式（529）是不成立的； 利用线性环节的脉冲传递函数只能得出在采样时刻上的信息。为了强调这一点，往往在环节的输出端画上一个假想的同步理想开关，如图58（b）所示。实际上，线性环节的输出仍然是一个连续信号c(t)。
脉冲传递函数的求法，连续系统的脉冲传递函数G(z)，可以通过其传递函数G(s)求取。具体步骤如下： 
(1) 对连续系统或元件的传递函数G(s)取拉普拉斯逆变换，求得脉冲响应g(t)为

g(t)=L－1［G(s)］


(2) 对g(t)进行Z变换，则得到系统或元件的脉冲传递函数G(z)。
脉冲传递函数还可由连续系统的传递函数，经部分分式法，通过查Z变换表求得。
MATLAB提供了如表51所示传递函数、脉冲传递函数、状态空间模型的仿真函数。


表51传递函数、脉冲传递函数、状态空间模型的仿真函数


函数
使 用 说 明
SYS=TF(NUM,DEN)
返回变量SYS为连续系统传递函数模型
SYS=TF(NUM,DEN,TS)
返回变量SYS为离散系统传递函数模型。TS为采样周期，当TS＝－1或者TS＝［］时，表示系统采样周期未定义
S=TF('s')
定义拉普拉斯变换算子 (拉普拉斯 variable)，以原形式输入传递函数
Z=TF('z',TS) 
定义Z变换算子及采样时间TS，以原形式输入传递函数
PRINTSYS(NUM,DEN,'s')
将系统传递函数以分式的形式打印出来，'s'表示传递函数变量续表


函数
使 用 说 明
PRINTSYS(NUM,DEN,'z')
将系统传递函数以分式的形式打印出来，'z'表示传递函数变量
GET(sys)
可获得传递函数模型对象sys的所有信息
SET(sys,'Property',Value,…)
为系统不同属性设定值
［NUM,DEN］=TFDATA(SYS,'v')
以行向量的形式返回传递函数分子分母多项式
C=CONV(A, B)
多项式A, B以系数行向量表示，进行相乘。结果C仍以系数行向量表示
sys=zpk(z,p,k)
得到连续系统的零极点增益模型
sys=zpk(z,p,k,Ts)
得到连续系统的零极点增益模型，采样时间为Ts 
s=zpk('s')
得到拉普拉斯算子，按原格式输入系统，得到系统zpk模型
z=zpk('z',Ts)
得到Z变换算子和采样时间Ts，按原格式输入系统，得到系统zpk模型
［Z,P,K］=ZPKDATA(SYS,'v')
得到系统的零极点和增益，参数'v'表示以向量形式表示
［p,z］ =pzmap(sys)
返回系统零极点
pzmap(sys)
得到系统零极点分布图
sys=ss(A,B,C,D)
由A，B，C，D矩阵直接得到连续系统状态空间模型
sys=ss(A,B,C,D,Ts)
由A，B，C，D矩阵和采样时间Ts直接得到离散系统状态空间模型
［A,B,C,D］ =ssdata(sys)
得到连续系统参数
［A,B,C,D,Ts］ =ssdata(sys)
得到离散系统参数

MATLAB提供的传递函数、脉冲传递函数、状态空间模型之间的转换函数见表52。


表52传递函数、脉冲传递函数、状态空间模型之间的转换函数


函数名
使 用 说 明
tfsys=tf(sys)
将其他类型的模型转换为多项式传递函数模型
zsys=zpk(sys)
将其他类型的模型转换为zpk模型
sys_ss=ss(sys)
将其他类型的模型转换为ss模型
［A,B,C,D］=tf2ss(num,den)
tf模型参数转换为ss模型参数
［num,den］=ss2tf(A,B,C,D,iu)
ss模型参数转换为tf模型参数，iu表示对应第i路传递函数
［z,p,k］=tf2zp(num,den)
tf模型参数转换为zpk模型参数
［num,den］=zp2tf(z,p,k)
zpk模型参数转换为tf模型参数
［A,B,C,D］=zp2ss(z,p,k)
zpk模型参数转换为ss模型参数
［z,p,k］=ss2zp(A,B,C,D,i)
ss模型参数转换为zpk模型参数，iu表示对应第i路传递函数
sys_min=minreal(sys)
对传递函数sys进行约分后，输出最小实现系统
sysd=c2d(sysc,Ts)
将连续系统转换为采样周期为Ts的离散系统
sysd=c2d(sysc,Ts, 'method')
指定连续系统的离散化方法
［sysd,G］ =c2d(sysc,Ts, 'method')
对于SS模型，求得初始条件的转换阵G
［Ad,Bd,Cd,Dd］ =c2dm(A,B,C,D,Ts,'method')
连续SS模型的离散化续表


函数名
使 用 说 明
sysc=d2c(sysd)
将离散系统转换为连续系统
sysc=d2c(sysd,method)
指定离散系统的连续化方法method
［Ac,Bc,Cc,Dc］ =d2cm(A,B,C,D,Ts,'method')
用于离散SS模型的连续化
sysd1 =d2d(sysd,Ts)
改变采样周期，生成新的离散系统

MATLAB提供的模型特征分析的相关函数见表53。


表53模型特征分析的相关函数


函数名
使 用 说 明
class
返回模型的类型名称，'tf'、'zpk'、'ss' 或者'frd'
hasdelay
如果LTI模型有任何类型的滞后时间，则返回1
isa
判断LTI模型的输入参量是否为指定的类型
isct
判断模型是否为连续系统模型
isdt
判断模型是否为离散系统模型
isempty
判断LTI模型是否为空
isproper
判断模型是否为正则系统(传递函数分母阶次大于或等于分子的阶次)
issiso
判断模型是否为MIMO系统
ndims
返回LTI数组的维数
reshape
改变LTI数组的形状
size
返回输入输出状态的维数
cover
计算输出的协方差和状态的协方差
damp
求取系统特征根的无阻尼自振频率和阻尼比
dcgain
返回系统的低频增益
dsort
离散系统的极点按幅值排序
esort
连续系统的极点按实部排序
norm
LTI模型的范数
pole，eig
求系统的极点
pzmap
求取系统的零极点分布
zero
LIT系统的传输零点

【例56】已知传递函数模型G（s）=10(2s+1)s2(s2+7s+13)
，将其输入MATLAB工作空间中。
方法1，程序如下： 

>>num=conv(10,［2,1］);

>>den=conv(［1 0 0］,［1 7 13］);

>>G=tf(num,den)

运行结果如下： 

G =

20 s + 10

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

s^4 + 7 s^3 + 13 s^2

 Continuous-time transfer function

方法2，程序如下： 

>>s=tf('s');

>>G=10*(2*s+1)/s^2/(s^2+7*s+13) 

运行结果如下： 

G =

20 s + 10

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

s^4 + 7 s^3 + 13 s^2

 Continuous-time transfer function

【例57】设置传递函数模型G（s）=10(2s+1)s2(s2+7s+13)
，时间延迟常数为τ=1．2，即系统模型为G(s)e-1.2s
。
可以通过在已有MATLAB模型基础上设置时间延迟常数来完成模型的输入，程序如下:

>> s=tf('s');

>> G=10*(2*s+1)/s^2/(s^2+7*s+13)； 

>> set(G,'ioDelay',1.2) %或 G.ioDelay=1.2

>> G

运行结果如下： 

G =

20 s + 10

exp(-1.2*s) * --------------------

s^4 + 7 s^3 + 13 s^2

 Continuous-time transfer function

【例58】已知一系统的传递函数G(s)=7s2+2s+84s3+12s2+4s+2
，求取其零极点向量和增益值，并得到系统的零极点增益模型。
实现程序如下： 

>> Gtf=tf(［7 2 8］,［4 12 4 2］)； 

>> ［z,p,k］=zpkdata(Gtf,'v')； 

>> Gzpk=zpk(z,p,k)

>> ［p1,z1］ = pzmap(Gtf) 

运行结果如下： 

Gzpk =

1.75 (s^2 + 0.2857s + 1.143)

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

(s+2.698) (s^2 + 0.302s + 0.1853)

Continuous-time zero/pole/gain model

p1 =

-2.6980 + 0.0000i

-0.1510 + 0.4031i

-0.1510 - 0.4031i

z1 =

-0.1429 + 1.0595i

-0.1429 - 1.0595i

【例59】已知系统的零极点模型G（s）=5(s+2)(s+4)(s+1)(s+3)
，求其TF模型及状态空间模型。
实现程序如下： 

>>z=［-2 -4］';

>>p=［-1 -3］';

>>k=5;

>>Gzpk=zpk(z,p,k);

>>［a,b,c,d］=zp2ss(z,p,k)

>>［num,den］=zp2tf(z,p,k)

>> Gtf=tf(Gzpk)

运行结果如下： 

a =

-4.0000-1.7321

1.73210

b =

1

0

c =

10.000014.4338

d =

5

num =

53040

den =

143

Gtf =

5 s^2 + 30 s + 40

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

s^2 + 4 s + 3

 Continuous-time transfer function

【例510】系统的被控对象传递函数为G（s）=10(s+2)(s+5)
，采样周期Ts＝0．01s，试将其进行离散化。
实现程序如下： 

>>num=10;

>>den=conv(［1,2］,［1,5］);

>>ts=0.01;

>>sysc=tf(num,den);

>>sysd=c2d(sysc,ts) 

运行结果如下： 

sysd = 

0.0004885 z + 0.0004772

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

z^2 - 1.931 z + 0.9324

 Sample time: 0.01 seconds

Discrete-time transfer function

2． 结构图模型
传递函数是由代数方程组通过消去系统中间变量得到的，如果系统结构复杂，方程组数目较多，那么消去中间变量就比较麻烦，并且中间变量的传递过程在系统输入与输出关系得不到反映，因此，结构图作为一种数学模型，在控制理论中得到了广泛的应用。结构图是将方块图与传递函数结合起来的一种将控制系统图形化了的数学模型。如果把组成系统的各个环节用方块表示，在方块内标出表征此环节输入输出关系的传递函数，并将环节的输入量、输出量改用拉普拉斯变换来表示，这种图形成为动态结构图，简称结构图。如果按照信号的传递方向将各环节的结构图依次连接起来，形成一个整体，这就是系统结构图。结构图不但能清楚表明系统的组成和信号的传递方向，而且能清楚地表示出系统信号传递过程中的数学关系。
建立系统结构图的步骤如下： 
(1) 首先应分别列写系统各元件的微分方程，在建立微分方程时，注意分清输入量和输出量，同时应考虑相邻元件之间是否存在负载效应。
(2) 设初始条件为零时，将各元件的微分方程进行取拉普拉斯变换，并作出各元件的结构图。
(3) 将系统的输入量放在最左边，输出量放在最右边，按照各元件的信号流向，用信号线依次将各元件的结构图连接起来，便构成系统的结构图。
MATLAB提供结构图模型的化简涉及的函数如表54所示


表54结构图模型的相关操作函数


函数名
使 用 说 明
sys=parallel(sys1,sys2)
并联两个系统，等效于sys=sys1+sys2
sys=parallel(sys1,sys2,inp1,inp2,out1,out2)
对MIMO系统，表示sys1的输入inp1与sys2的输入inp2相连，sys1输出out1与sys2输出out2相连
sys=series(sys1,sys2)
串联两个系统，等效于sys=sys2 * sys1
sys=feedback(sys1,sys2)
两系统负反馈连接，默认格式
sys=feedback(sys1,sys2,sign)
sign=-1表示负反馈，sign=1表示正反馈。等效于sys=sys1/(1±sys1*sys2) 
sys=feedback(sys1,sys2,feedin,feedout,sign)
对MIMO系统，部分反馈连接。sys1的指定输出feedout连接到sys2的输入，而sys2的输出连接到sys1的指定输入feedin，最终实现的反馈系统与sys1具有相同的输入、输入端。sign标识正负反馈
sys=append(sys1,sys2,…,sysN)
将子系统sys1，sys2，…，sysN的所有输入作为系统的输入，所有输出作为系统输出，且各子系统间没有信号连接，从而扩展为一个系统
sysc=connect(sys, Q, inputs, outputs)
将多个子系统按照一定的连接方式构成一个系统。sys是待连接的子系统被append函数扩展后的系统。Q矩阵声明了子系统的连接方式。Q矩阵的行向量声明了sys输入信号的连接方式，每个行向量的第1个元素为sys系统的输入端口号，其他元素为与该输入信号相连接的sys端口号。inputs声明了整个系统的输入信号是由sys系统的哪些输入端口号构成。outputs声明了整个系统的输出信号是由sys系统的哪些输出端口号构成

【例511】化简如图59所示的系统，求系统的传递函数。


图59例511系统图


实现程序如下： 

>>clear

>>G1=tf(1,［1 1］);

>>G2=tf(1,［3 4 1］);

>>Gp=G1+G2;

>>G3=tf(1,［1 0］);

>>Gs=series(G3,Gp);

>>Gc=Gs/(1+Gs)

>>Gc1=minreal(Gc) 

运行结果如下： 

Gc =

9 s^6 + 36 s^5 + 56 s^4 + 42 s^3 + 15 s^2 + 2 s

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

9 s^8 + 42 s^7 + 88 s^6 + 112 s^5 + 95 s^4 + 52 s^3 + 16 s^2 + 2 s

 Continuous-time transfer function.

Gc1 =

s^4 + 3.667 s^3 + 5 s^2 + 3 s + 0.6667

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

s^6 + 4.333 s^5 + 8.333 s^4 + 9.667 s^3 + 7.333 s^2 + 3.333 s + 0.6667

 Continuous-time transfer function.

3． 信号流图模型
系统的结构图可以用来直接绘制信号流图，以图510为例简要说明信号流图的相关概念。


图510典型闭环控制系统的信号流图


(1) 源节点(输入节点)： 只有输出支路而没有输入支路的节点，称为源节点。它一般表示系统的输入变量，亦称输入节点，如图510中的节点R和N。
(2) 阱节点(输出节点)： 只有输入支路而没有输出支路的节点，称为阱节点。它一般表示系统的输出变量，亦称输出节点，如图510中的节点C。
(3) 混合节点： 既有输入支路又有输出支路的节点，称为混合节点，如图510中的节点E、Q、O。
(4) 通路： 沿着支路箭头的方向顺序穿过各相连支路的路径，如图510中的R、E、Q、O、C等。
(5) 前向通路： 从源节点出发并且终止于阱节点，与其他节点相交不多于一次的通路称为前向通路，如图510中的N、Q、O、C等。
(6) 回路： 起点和终点在同一个节点，并且与其他节点相交不多于一次的闭合路径称为回路，如图510中的E、Q、O、E。
(7) 前向通路传输(增益)： 前向通路中各支路传输(增益)的乘积。
(8) 回路传输(增益)： 回路中各支路传输(增益)的乘积。
(9) 不接触回路： 信号流图中，没有任何公共节点的回路，称为不接触回路或互不接触回路。
用信号流图代替系统的结构图，其优点在于不简化信号流图的情况下，利用梅逊公式直接求得源节点和阱节点之间的总增益。对于动态系统来说，这个总增益就是系统相应的输入和输出间的传递函数。
计算任意输入节点和输出节点之间传递函数G(s)的梅逊公式为

G(s)=1Δ·∑nk=1PkΔk(530)

式中，Δ——特征式，其计算公式为

Δ=1－∑La+∑LbLc－∑LdLeLf+…


n——从输入节点到输出节点间前向通路的条数； 
Pk——从输入节点到输出接点间第k条前向通路的总增益； 
∑La——所有不同回路的回路增益之和； 
∑LbLc——所有两两互不接触回路的回路增益乘积之和； 
∑LdLeLf——所有3个互不接触回路的回路增益乘积之和。
Δk——第k条前向通路的余子式，即把与该通路相接触的回路增益置为0后，特征式Δ所余下的部分。
【例512】根据典型闭环控制系统的信号流图，求传递函数C(s)/R(s)。
实现程序如下： 

>>syms g1 g2 H

>>syms R x1 x2 x3 x4 N

>> RaC=solve('x1=R-x3*H','x2=x1*g1','x3=x2*g2','x4=x3','x1,x2,x3,x4'); 

>>pretty(simplify(RaC.x4/R))

运行结果如下： 

g1 g2

-----------

H g1 g2 + 1

4． 频率特性模型
设线性定常系统的传递函数为G(s)，其对正弦输入信号的稳态响应仍然是与输入信号同频率的正弦信号，输出信号的振幅是输入信号的|G(jω)|倍，输出信号相对输入信号的相移为φ=∠G(jω)，输出信号的振幅及相移都是角频率ω的函数。
G(jω)=|G(jω)|ej∠G(jω)称为系统的频率特性，它表明了正弦信号作用下，系统的稳态输出与输入信号的关系。其中，|G(jω)|为幅频特性，它反映了系统在不同频率的正弦信号作用下，稳态输出的幅值与输入信号幅值之比。∠G(jω)=arctanImG(jω)ReG(jω)称为相频特性，它反映了系统在不同频率的正弦信号作用下，输出信号相对输入信号的相移。系统的幅频特性和相频特性统称为系统的频率特性。
比较系统的频率特性和传递函数可知，频率特性与传递函数有如下关系：

G(jω)=G(s)|s=jω(531)

一般地，若系统具有以下传递函数： 

G(s)=bmsm+bm－1sm－1+…+b1s+b0ansn+an－1sn－1+…+a1s+a0(532)

则系统频率特性可写为

G(jω)=bm(jω)m+bm－1(jω)m－1+…+b1(jω)+b0an(jω)n+an－1(jω)n－1+…+a1(jω)+a0(533)


由式(532)可推导出线性定常系统的频率特性。对于稳定的系统，可以由实验的方法确定系统的频率特性，即在系统的输入端作用不同频率的正弦信号，在输出端测得相应的稳态输出的幅值和相角，根据幅值比和相位差，就可得到系统的频率特性。对于不稳定的系统，则不能由实验的方法得到系统的频率特性，这是由于系统传递函数中不稳定极点会产生发散或振荡的分量，随时间推移，其瞬态分量不会消失，所以不稳定系统的频率特性是观察不到的。
频率特性的曲线表示方法主要有3种形式，即对数坐标图、极坐标图、对数幅相图。
对数坐标图又称伯德(Bode)图，包括对数幅频和对数相频两条曲线。在实际应用中，经常采用这种曲线来表示系统的频率特性。对数幅频特性曲线的横坐标是频率ω，按对数分度，单位是rad/s。纵坐标表示对数幅频特性的函数值，采用线性分度，单位是dB。对数幅频特性用L(ω)表示，定义为： L(ω)=20log|G(jω)|，对数相频特性曲线的横坐标也是频率ω，按对数分度，单位是rad/s。纵坐标表示相频特性的函数值，记作φ(ω)，单位是度。
极坐标图又称幅相频率特性曲线、奈奎斯特曲线。其特点是以角频率ω作自变量，把幅频特性和相频特性用一条曲线同时表示在复平面上。
对数幅相图又称Nichols曲线，是将对数幅频特性和对数相频特性两张图，在角频率ω为参变量的情况下合成一张图。即以相位φ(ω)为横坐标，以20logA(ω)为纵坐标，以ω为参变量的一种图示法。
MATLAB提供的频率特性相关的函数如表55所示，注意bode、nyquist、nichols函数前加上的d，可得到绘制离散系统的对应功能图的函数，一般还要指定采样周期。


表55频率特性相关函数的说明


函数名
使 用 说 明
bode(G)
绘制系统伯德图
bode(G,w)
绘制指定频率范围的系统伯德图续表


函数名
使 用 说 明
bode(G1,'r--',G2,'gx',…)
同时绘制多系统伯德图。图形属性参数可选
［mag,phase,w］=bode(G)
返回系统伯德图相应的幅值、相位和频率向量，
可用magdb=20*log10(mag)将幅值转换为分贝值
［mag,phase］=bode(G,w)
返回系统伯德图与指定w相应的幅值、相位
nyquist(sys)
绘制系统奈奎斯特图。系统自动选取频率范围
nyquist(sys,w)
绘制系统奈奎斯特图。由用户指定选取频率范围
nyquist(G1,'r--',G2,'gx',…)
同时绘制多系统带图形参数奈奎斯特图
［re,im,w］=nyquist(sys)
返回系统奈奎斯特图相应的实部、虚部和频率向量
［re,im］=nyquist(sys,w)
返回系统奈奎斯特图与指定w相应的实部、虚部
nichols(G)
绘制系统Nichols图
nichols(G,w)
绘制指定频率范围的系统Nichols图
nichols(G1,'r--',G2,'gx',…)
同时绘制多系统带图形参数Nichols图
［mag,phase,w］=nichols(G)
返回系统Nichols图相应的幅值、相位和频率向量
［mag,phase］=nichols(G,w)
返回系统Nichols图与指定w相应的幅值、相位
sigma(sys)
绘制系统sys的奇异值曲线
evalfr(sys,x)
计算系统sys的单个复频率点x的频率响应
freqresp(sys,w)
计算系统sys在给定频率区间w的频率响应
ngrid
在Nichols曲线图上绘制等M圆和等N圆
ngrid('new')
绘制网格前清除原图，然后设置hold on

【例513】二阶振荡环节的传递函数为G(s)=ω2ns2+2ζωns+ω2n，绘制ζ取不同值时的伯德图，ωn=5，ζ取0．1∶0．2∶1．0。
实现程序如下：

>>wn=5;

>>kosi=［0.1:0.2:1.0］;

>>w=logspace(-1,1,100);

>>figure(1)

>>num=［wn.^2］;

>>for kos=kosi 

>>den=［1 2*kos*wn wn.^2］;

>>［mag,pha,w1］=bode(num,den,w);

>>subplot(2,1,1);hold on

>>semilogx(w1,mag);

>>subplot(2,1,2);hold on

>>semilogx(w1,pha);

>>end

>>subplot(2,1,1);grid on

>>title('Bode Plot');

>>xlabel('Frequency(rad/sec)');

>>ylabel('Gain dB');

>>subplot(2,1,2);grid on

>>xlabel('Frequency(rad/sec)');

>>ylabel('Phase deg');

>>hold off

运行结果如图511所示。


图511例513程序运行结果图







【例514】系统的传递函数为G(s)=30(10s+1)(3s+1)(0．5s+1)(0．03s+1)(s2+s+1)，采用MATLAB绘制系统的渐进幅频特性。
先编写asympbode．m函数用于计算绘制渐进幅频特性的数据。

function ［wpos,ypos］=asympbode(G,w)

G1=zpk(G);

wpos=［］;

pos1=［］;

if nargin==1,w=freqint2(G);

end

zer=G1.z{1}; pol=G1.p{1};

gain=G1.k;

for i=1:length(zer);

ifisreal(zer(i))

wpos=［wpos,abs(zer(i))］;

pos1=［pos1,20］;

else

if imag(zer(i))>0

wpos=［wpos,abs(zer(i))］;

pos1=［pos1,40］;

end

end

end

for i=1:length(pol);

if isreal(pol(i))

wpos=［wpos,abs(pol(i))］;

pos1=［pos1,-20］;

else

if imag(pol(i))>0

wpos=［wpos,abs(pol(i))］;

pos1=［pos1,-40］;

end

end

end

wpos=［wpos w(1) w(length(w))］;

pos1=［pos1,0,0］;

［wpos,ii］=sort(wpos);

pos1=pos1(ii);

ii=find(abs(wpos)<eps);

kslp=0;

w_start=1000*eps;

if length(ii)>0

kslp=sum(pos1(ii));

ii=(ii(length(ii))+1):length(wpos);

wpos=wpos(ii);

pos1=pos1(ii);

end

while 1

［ypos1,pp］=bode(G,w_start);

if isinf(ypos1),w_start=w_start*10;

else break;

end

end

wpos=［w_start wpos］;

ypos(1)=20*log10(ypos1);

pos1=［kslp pos1］;

for i=2:length(wpos)

kslp=sum(pos1(1:i-1));

ypos(i)=ypos(i-1)+kslp*log10(wpos(i)/wpos(i-1));

end

ii=find(wpos>=w(1)&wpos<=w(length(w)));

wpos=wpos(ii);

ypos=ypos(ii);

然后在命令窗口中输入： 

>> s=tf('s');

>>G1=30*(10*s+1)/((3*s+1)*(0.5*s+1)*(0.03*s+1)*(s*s+s+1));

>> w=10e-3:0.1:1000;

>>［x1,y1］=asympbode (G1,w);

>>semilogx(x1,y1)

>> grid;

运行结果如图512所示。


图512例514程序运行结果图


5．3神经网络建模及示例
神经网络是由大量人工神经元(处理单元)广泛互联而成的网络，它是在现代神经生物学和认识科学对人类信息处理研究的基础上提出来的，具有很强的自适应性和学习能力、非线性映射能力、鲁棒性和容错能力。充分地将这些神经网络特性应用于控制领域，可使控制系统的智能化向前迈进一大步。随着被控系统越来越复杂，人们对控制系统的要求越来越高，特别是要求控制系统能适应不确定的、时变的对象与环境。传统的基于精确模型的控制方法难以适应要求，现在关于控制的概念也已更加广泛，要求包括一些决策、规划以及学习功能。神经网络由于具有上述优点而越来越受到人们的重视。本书介绍最为常用的BP(Back Propagation)网络。
人工神经元是神经网络的基本元素，其原理可以用图513表示。


图513人工神经元结构示意图

人工神经元是对生物神经元的一种模拟与简化，它是神经网络的基本处理单元。如图所示为一种简化的人工神经元结构。它是一个多输入、单输出的非线性元件。
其输入、输出关系为

Ii=∑nj=1wjixj－θi

yi=f(Ii)

其中，xj(j=1,2,…,n)是从其他神经元传来的输入信号； wji表示神经元i到神经元j的连接权值； θi为阈值； f(·)称为激发函数或作用函数。方便起见，常把－θi看成是恒等于1的输入x0的权值，则上式可写成

Ii=∑nj=0wjixj

其中，w0i=－θi，x0=1。
输出激发函数f(·)又称为变换函数，它决定神经元(节点)的输出。该输出为0或1，取决于其他输入之和是大于还是小于内部阈值θi。f(·)一般具有非线性特征。常用的激活函数有线性函数、斜坡函数、阈值函数，这3个激活函数都属于线性函数，还有两个常用的非线性激活函数。分别是S形函数和双极S形函数。双极S形函数与S形函数的主要区别在于函数的值域，双极S形函数值域是(-1,1)，而S形函数值域是(0,1)。由于S形函数与双极S形函数都是可导的，因此适合用在BP神经网络中。根据网络中神经元的互联方式，常见网络结构主要可以分为如下３类。
前馈神经网络： 前馈网络也称前向网络。这种网络只在训练过程会有反馈信号，而在分类过程中数据只能向前传送，直到到达输出层，层间没有向后的反馈信号，因此被称为前馈神经网络。感知机与BP神经网络就属于前馈网络。
反馈神经网络： 反馈型神经网络是一种从输出到输入具有反馈连接的神经网络，其结构比前馈网络要复杂得多。典型的反馈型神经网络有Elman网络和Hopfield网络。
自组织网络： 自组织神经网络是一种无导师学习网络。它通过自动寻找样本中的内在规律和本质属性，自组织、自适应地改变网络参数与结构。
误差反向传播神经网络，简称BP网络，是一种单向传播的多层前向网络。在模式识别、图像处理、系统辨识、函数拟合、优化计算、最优预测和自适应控制等领域有着较为广泛的应用。图514是BP网络的示意图。


图514BP网络示意图


误差反向传播的BP算法简称BP算法，其基本思想是最小二乘算法。采用梯度搜索技术，以使网络的实际输出值与期望输出值的误差均方值为最小。BP算法的学习过程由正向传播和反向传播组成。在正向传播过程中，输入信息从输入层经隐含层逐层处理，并传向输出层，每层神经元(节点)的状态只影响下一层神经元的状态。如果在输出层不能得到期望的输出，则转入反向传播，将误差信号沿原来的连接通路返回，通过修改各层神经元的权值，使误差信号最小。
设BP网络的结构如图所示，有M个输入节点，输入层节点的输出等于其输入。输出层有L个输出节点，网络的隐含层有q个节点，w ij是输入层和隐含层节点之间的连接权值。w jk是隐含层和输出层节点之间的连接权值，隐含层和输出层节点的输入是前一层节点的输出的加权和，每个节点的激励程度由它的激发函数来决定。
1． BP网络的前馈计算
在训练该网络的学习阶段，设有N个训练样本，先假定用其中的某一固定样本中的输入/输出模式对网络进行训练。若网络输出与期望输出值不一致，则将其误差信号从输出端反向传播，并在传播过程中对加权系数不断修正，使在输出层节点上得到的输出结果尽可能接近期望输出值。对样本完成网络加权系数的调整后，再送入另一样本模式对，进行类似的学习，直到完成所有样本的训练学习为止。
2． BP网络权值的调整规则
设每一样本p的输入输出模式对的二次型误差函数定义为

Ep=12∑Lk=1(dpk－Opk)2


系统的平均误差代价函数为

E=12∑Pp=1∑Lk=1(dpk－Opk)2=12∑Pp=1Ep

式中，P为样本模式对数，L为网络输出节点数。
下面介绍基于一阶梯度法的优化方法，即最速下降法。简便起见，略去下标p，有

E=12∑Lk=1(dk－Ok)2


权系数应按E函数梯度变化的反方向进行调整，使网络的输出接近期望的输出。
(1) 输出层权系数的调整。
权系数的修正公式为

Δwjk=－ηEwjk

式中，η为学习速率，η>0； 

Ewjk=Enetknetkwjk


定义反传误差信号δk为

δk=－Enetk=EOkOknetk

式中： 

EOk=－(dk－Ok)
Oknetk=netkf(netk)=f′(netk)

因此

δk=(dk－Ok)f′(netk)=Ok(1－Ok)(dk－Ok)
netkwjk=wjk∑qj=1wjkOj=Oj

由此可得输出层的任意神经元权系数的修正公式为

Δwjk=η(dk－Ok)f′(netk)Oj=ηδkOj

=ηOk(1－Ok)(dk－Ok)Oj


(2) 隐含层节点权系数的调整。
计算权系数的变化量为

Δwij=－ηEwij=－ηEnetjnetjwij=－ηEnetjOi

=η－EOjOjnetjOi=η－EOjf′(netj)Oi=ηδjOi

式中，E/Oj不能直接计算，需通过其他间接量进行计算，即

－EOj=－∑Lk=1EnetknetkOj=∑Lk=1－EnetkOj∑qj=1wjkOj

=∑Lk=1－Enetkwjk=∑Lk=1δkwjk

显然有

δj=f′(netj)∑Lk=1δkwjk


将样本标记p记入公式后，对于输出节点k，有

Δpwjk=ηf′(netpk)(dpk－Opk)Opj=ηOpk(1－Opk)(dpk－Opk)Opj

对于隐含节点j，有

Δpwij=ηf′(netpj)∑Lk=1δpkwjkOpi=ηOpj(1－Opj)∑Lk=1δpkwjkOpi

式中，O pk是输出节点k的输出； O pj是隐含节点j的输出； O pi是输入节点i的输出。
从上面推导的结果可得网络连接权值调整式

wij(t+1)=wij(t)+ηδiOi+α［wij(t)－wij(t－1)］



图515BP学习算法流程图

式中，t+1表示第t+1步； α为平滑因子，0<α<1。
3． BP学习算法的计算步骤
（1） 初始化： 置所有权值为较小的随机数。
（2） 提供训练集： 给定输入向量X=(x1,x2,…,xM)和期望的目标输出向量D=(d0,d1,…,dL)。
（3） 计算实际输出，计算隐含层、输出层各神经元输出。
（4） 计算目标值与实际输出的偏差E。
（5） 计算Δpwjk。
（6） 计算Δpwij。
（7） 返回步骤（2）重复计算，直到误差Ep满足要求为止。
BP学习算法流程图如图515所示。
4． 使用BP算法应注意的问题
（1） 隐含层节点的个数对于识别率的影响并不大，但是节点个数过多会增加运算量，使得训练较慢。
（2） 学习开始时，各隐含层连接权系数的初值应以设置较小的随机数较为适宜。
（3） 采用S型激发函数时，由于输出层各神经元的输出只能趋于1或0，不能达到1或0。在设置各训练样本时，期望的输出分量不能设置为1或0。
（4） 学习速率η的选择，在学习开始阶段，η选较大的值可以加快学习速度。学习接近优化区时，η值必须相当小，否则权系数将产生振荡而不收敛。
MATLAB提供了实现神经网络功能的大量函数和集成化GUI工具箱nntool，表56给出了神经网络相关的函数说明，供编写复杂神经网络程序时查用。


表56神经网络的相关函数说明


函数名
使 用 说 明
函数名
使 用 说 明
函数名
使 用 说 明
newp
创建感知器网络
newlind
设计线性层
newlin
创建线性层
newff
创建前馈BP网络
newcf
创建多层前馈BP网络
newfftd
创建前馈输入延迟BP网络
newrb
创建径向基网络
newrbe
创建严格的径向基网络
newgrnn
创建广义回归神经网络
newpnn
创建概率神经网络
newc
创建竞争层
newsom
创建自组织特征映射
sim
仿真神经网络
init
初始化神经网络
adapt
神经网络的自适应化
train
训练神经网络
dotprod
权函数的点积
ddotprod
权函数点积的导数
dist
Euclidean距离权函数
normprod
规范点积权函数
negdist
Negative距离权函数
mandist
Manhattan距离权函数
linkdist
Link距离权函数
netsum
网络输入函数的求和
dnetsum
网络输入函数求和的导数
hardlim
硬限幅传递函数
hardlims
对称硬限幅传递函数
purelin
线性传递函数
tansig
正切S形传递函数
logsig
对数S形传递函数
dpurelin
线性传递函数的导数
dtansig
正切S形传递函数的导数
dlogsig
对数S形传递函数的导数
compet
竞争传递函数
radbas
径向基传递函数
satlins
对称饱和线性传递函数
initlay
层与层之间的网络初始化函数
initwb
阈值与权值的初始化函数
initzero
零权/阈值的初始化函数
initnw
Nguyen_Widrow层的初始化函数
initcon
Conscience阈值的初始化函数
midpoint
中点权值初始化函数
mae
均值绝对误差性能分析函数
mse
均方差性能分析函数
msereg
均方差w/reg性能分析函数
dmse
均方差性能分析函数的导数
dmsereg
均方差w/reg性能分析函数的导数
learnp
感知器学习函数
learnpn
标准感知器学习函数
learnwh
Widrow_Hoff学习规则
learngd
BP学习规则
learngdm
带动量项的BP学习规则
learnk
Kohonen权学习函数
learncon
Conscience阈值学习函数
learnsom
自组织映射权学习函数
adaptwb
网络权与阈值的自适应函数
trainwb
网络权与阈值的训练函数
traingd
梯度下降的BP算法训练函数
traingdm
梯度下降w/动量的BP算法训练函数
traingda
梯度下降w/自适应lr的BP 算法训练函数续表


函数名
使 用 说 明
函数名
使 用 说 明
函数名
使 用 说 明
traingdx
梯度下降w/动量和自适应lr的BP算法训练函数
trainlm
Levenberg_Marquardt的BP算法训练函数
trainwbl
每个训练周期用一个权值向量或偏差向量的训练函数
maxlinlr
线性学习层的最大学习率
gridtop
网络层拓扑函数
randtop
随机层拓扑函数

如果要使用集成化GUI工具箱nntool来完成神经网络模型的建立，只需在MATLAB命令窗口中命令行输入nntool，就会运行神经网络工具箱的主界面，主界面由6个部分组成： 系统的输入数据、系统的期望输出、网络的计算输出、网络的误差、已经建立的神经网络以及数据的导入和网络模型的建立。下面通过一个例子来说明其使用方法。
【例515】采用神经网络工具箱nntool来逼近函数f（x）=0．08sincos（x）+e-x，x∈［1,20］。
第1步，在MATLAB命令行窗口中输入nntool后按回车键，打开GUI编辑窗口，如图516所示。


图516Network/Data Manager窗口


其各区域和按钮的含义如下： 
 Input Data区域——显示指定的输入数据。
 Output Data区域——显示网络对应的输出数据。
 Target Data区域——显示期望输出数据。
 Error Data区域——显示误差。
 Input Delay States区域——显示设置的输入延迟参数。
 Layer Delay States区域——显示层的延迟状态。
 Networks区域——显示设置网络的类型。
 Help按钮——Network/Data Manager窗口有关区域和按钮的详细介绍。
 Close按钮——关闭Network/Data Manager窗口。
 Import按钮——将MATLAB工作空间或者文件中的数据和网络导入GUI工作空间。
 New按钮——创建新网络或者生成新的数据。
 Open按钮——打开选定的数据或网络，以便查看和编辑。
 Export按钮——将GUI工作空间的数据和网络导出到MATLAB工作空间或文件中。
 Delete按钮——删除所选择的数据或者网络。
第2步，导入数据，在主界面Neural Network/Data Manager窗口中单击Import按钮，弹出Import to Network/Data Manager 窗口，如图517所示。


图517Import to Network/Data Manager窗口


可以在左侧Source选项组中选择从MATLAB工作空间导入数据或者加载文件中的数据，这里选择从工作空间导入。首先在MATLAB工作空间中输入如下语句来模拟未知系统的输入向量x和目标输出向量y。

>> x=［1:0.05:20］;

>>y=0.08.*sin(cos(x))+exp(-x);

此时Import to Network/Data Manager窗口中的Select a Variable列表框显示出当前工作空间中的所有变量，从中选择输入向量x，并在Import As选项组中选中Input Data，单击Import按钮； 再选择目标输出向量y，并在Import As列表框中选中Target Data，单击Import按钮，就完成了输入输出数据的导入。单击Close按钮返回Neural Network/Data Manager主窗口。在主窗口中双击向量x，就可以查看变量的内容，同理，向量y也一样能查看到。
第3步，创建BP网络。在Neural Network/Data Manager主窗口，单击New按钮，在弹出的Create Network or Data窗口中选择Network选项，开始创建网络。根据要求在窗口中修改相应的参数： 在Name 文本框中输入所创建的网络名称，将该网络命名为BPexenet。在Network Type下拉列表框中选择需要创建的网络类型，此处为Feedforward backprop。在Input data和Target data的下拉列表框中选择相应的导入数据。Training function下拉列表框用于设置该网络训练函数的类型，这里选择TRAINLM。在Adaption learning function下拉列表框中设置网络的学习函数LEARNGDM，网络的性能函数Performance function取默认值MSE，网络的层数Number of layers输入2，表明所设计的网络有一层是隐含层。最后设置网络各层的传递函数和神经元的数目。在Properties for下拉列表框中选择Layer1，表示接下来设置隐含层的传递函数和神经元的数目。在Number of neurons 文本框中输入20，表示隐含层由20个神经元组成。在Transfer Function下拉列表框中选择传递函数类型为TANSIG，如图518所示。然后在Properties for下拉列表框中选择Layer2，表示接下来设置输出层属性，输出层的Number of neurons文本框中不能输入数字，网络自动默认为1，在Transfer Function下拉列表框中选择传递函数类型为PURELIN。单击View按钮可以看到所建立BP神经网络的结构，如图519所示。单击图518中的Create按钮，就完成了神经网络的设置。单击图518中的Close按钮返回Neural Network/Data Manager主窗口。



图518网络创建窗口




图519BPexenet的网络的结构


第4步，训练网络。在Neural Network/Data Manager窗口中选中生成的BP神经网络BPexenet，则Open、Delete按钮被激活。单击Open按钮，该窗口用于设置网络自适应、仿真、训练和重新初始化的参数，并执行相应的过程。单击Train选项卡后做相应的设置即可进行神经网络的训练： 在Training Info选项卡中，选择输入Inputs为P，期望输出Targets为T。在Training Parameters选项卡中，设置训练步数epochs为100，目标误差goal为0．001，其他参数采用默认值，如图520所示。


图520设置训练参数窗口


参数修改完毕后，单击Network： BPexenet窗口右下角的Train Network按钮，就开始了网络训练，训练完成后会有一个结果信息界面，如图521所示。


图521网络训练结果


单击Performance、Training State以及Regression可以分别显示网络训练误差随训练次数的变化情况，训练数据的梯度、均方误差以及验证集检验失败的次数随训练次数的变化情况和训练集、验证集、测试集和全集的回归系数，回归系数越接近1表示模型训练的性能越好，如图522所示。


图522网路训练性能、网络训练状态、网络训练的回归分析


第5步，网络仿真。通过网络仿真可以验证训练后的神经网络对模型的逼近效果。在Network： BPexenet窗口中，单击Simulate标签，在Inputs下拉菜单中选择输入量x，同时将网络仿真输出Outputs变为BPexenetOutSim，以区别于训练时的输出BPexenetoutputs。单击Simulate Network按钮，回到Neural Network/Data Manager窗口，就会看到在Outputs区域有一个新的变量——BPexenetOutSim。双击变量名BPexenetOutSim，同样会弹出一个新的窗口，显示出对应输入向量x的网络输出值，主界面上单击Export按钮，就能在Output Data中选择想要查看、保存的数据，将其导出到MATLAB的工作空间。在MATLAB工作空间中输入语句plot(x,y,x,BPexenet_outSim)来验证模型的逼近效果，如图523所示。



图523对比结果


习题
1． 已知系统时域模型为y′=y2－t－24(t+1),0≤t≤1

y(0)=2试求其数值解，并与精确解相比较，其精确解为y(t)=t+1+1。
2． 求解差分方程y(n)-0．8y(n-1)-0．3y(n-2)=0．1x(n)+0．2x(n-1)+0．1x(n-2)，初值为y(-1)=0，y(-2)=1，x(-1)=1，x(-1)=1。
3． 设某离散系统的脉冲传递函数为G(z)=0．3z3+0．5z2+0．3z+0．8z4+3．2z3+3．9z2+2．2z+0．4，采样周期为T=0．01s，将其输入MATLAB的工作环境中，并且绘制零、极点分布图，并将该离散系统脉冲传递函数模型转换成状态空间表达式。
4． 已知线性定常系统的状态空间模型为

x·1
x·2=01
－1－2x1
x2+10
11u1
u2
y1
y2
y3=31
12
－2－1x1
x2+10
00
01u1
u2

试采用MATLAB求出其传递函数模型。
5． 典型反馈控制系统结构如图524所示，其中G(s)=0．5s2+0．8s4+3．2s3+3．9s2+2．2s+0．4，Gc(s)=s－1s+1，H(s)=10．1s+1，试采用MATLAB求出其传递函数。


图524习题5的系统结构


6． 采用神经网络工具箱nntool来逼近函数f(x)=30sine-xcos(x)+e-x，x∈［1,20］。