控制系统数字仿真习题答案.doc

第一章 习题
1-1什么是仿真？它所遵循的基本原则是什么？
答：仿真是建立在控制理论，相似理论，信息处理技术和计算技术等理论基础之上的，以计算机和其他专用物理效应设备为工具，利用系统模型对真实或假想的系统进行试验，并借助专家经验知识，统计数据和信息资料对试验结果进行分析和研究，进而做出决策的一门综合性的试验性科学。
它所遵循的基本原则是相似原理。
1-2在系统分析与设计中仿真法与解析法有何区别？各有什么特点？
答：解析法就是运用已掌握的理论知识对控制系统进行理论上的分析，计算。它是一种纯物理意义上的实验分析方法，在对系统的认识过程中具有普遍意义。由于受到理论的不完善性以及对事物认识的不全面性等因素的影响，其应用往往有很大局限性。
仿真法基于相似原理，是在模型上所进行的系统性能分析与研究的实验方法。
1-3数字仿真包括那几个要素？其关系如何？
答: 通常情况下，数字仿真实验包括三个基本要素，即实际系统，数学模型与计算机。由图可见，将实际系统抽象为数学模型，称之为一次模型化，它还涉及到系统辨识技术问题，统称为建模问题；将数学模型转化为可在计算机上运行的仿真模型，称之为二次模型化，这涉及到仿真技术问题，统称为仿真实验。
1-4为什么说模拟仿真较数字仿真精度低？其优点如何？。
答：由于受到电路元件精度的制约和容易受到外界的干扰，模拟仿真较数字仿真精度低
但模拟仿真具有如下优点：
（1） 描述连续的物理系统的动态过程比较自然和逼真。
（2） 仿真速度极快，失真小，结果可信度高。
（3） 能快速求解微分方程。模拟计算机运行时各运算器是并行工作的，模拟机的解题速度与原系统的复杂程度无关。
（4） 可以灵活设置仿真试验的时间标尺，既可以进行实时仿真，也可以进行非实时仿真。
（5） 易于和实物相连。
1-5什么是CAD技术？控制系统CAD可解决那些问题？
答：CAD技术，即计算机辅助设计（Computer Aided Design）,是将计算机高速而精确的计算能力，大容量存储和数据的能力与设计者的综合分析，逻辑判断以及创造性思维结合起来，用以快速设计进程，缩短设计周期，提高设计质量的技术。
控制系统CAD可以解决以频域法为主要内容的经典控制理论和以时域法为主要内容的现代控制理论。此外，自适应控制，自校正控制以及最优控制等现代控制测略都可利用CAD技术实现有效的分析与设计。
1-6什么是虚拟现实技术？它与仿真技术的关系如何？
答：虚拟现实技术是一种综合了计算机图形技术，多媒体技术，传感器技术，显示技术以及仿真技术等多种学科而发展起来的高新技术。
1-7什么是离散系统？什么是离散事件系统？如何用数学的方法描述它们？
答：本书所讲的“离散系统”指的是离散时间系统，即系统中状态变量的变化仅发生在一组离散时刻上的系统。它一般采用差分方程，离散状态方程和脉冲传递函数来描述。
离散事件系统是系统中状态变量的改变是由离散时刻上所发生的事件所驱动的系统。这种系统的输入输出是随机发生的，一般采用概率模型来描述。
1-8如图1-16所示某卫星姿态控制仿真实验系统，试说明：
（1） 若按模型分类，该系统属于那一类仿真系统？
（2） 图中“混合计算机”部分在系统中起什么作用？
（3） 与数字仿真相比该系统有什么优缺点？
答：（1）按模型分类，该系统属于物理仿真系统。
（2）混合计算机集中了模拟仿真和数字仿真的优点，它既可以与实物连接进行实时仿真，计算一些复杂函数，又可以对控制系统进行反复迭代计算。其数字部分用来模拟系统中的控制器，而模拟部分用于模拟控制对象。
（4） 与数字仿真相比，物理仿真总是有实物介入，效果逼真，精度高，具有实时性与在线性的特点，但其构成复杂，造价较高，耗时过长，通用性不强。

第二章习题
2-1 思考题：
（1）数学模型的微分方程，状态方程，传递函数，零极点增益和部分分式五种形式，各有什么特点？
（2）数学模型各种形式之间为什么要互相转换？
（3）控制系统建模的基本方法有哪些？他们的区别和特点是什么？
（4）控制系统计算机仿真中的“实现问题”是什么含意？
（5）数值积分法的选用应遵循哪几条原则？
答：（1）微分方程是直接描述系统输入和输出量之间的制约关系，是连续控制系统其他数学模型表达式的基础。状态方程能够反映系统内部各状态之间的相互关系，适用于多输入多输出系统。传递函数是零极点形式和部分分式形式的基础。零极点增益形式可用于分析系统的稳定性和快速性。利用部分分式形式可直接分析系统的动态过程。
（2）不同的控制系统的分析和设计方法，只适用于特定的数学模型形式。
（3）控制系统的建模方法大体有三种：机理模型法，统计模型法和混合模型法。机理模型法就是对已知结构，参数的物理系统运用相应的物理定律或定理，经过合理的分析简化建立起来的各物理量间的关系。该方法需要对系统的内部结构和特性完全的了解，精度高。统计模型法是采用归纳的方法，根据系统实测的数据，运用统计规律和系统辨识等理论建立的系统模型。该方法建立的数学模型受数据量不充分，数据精度不一致，数据处理方法的不完善，很难在精度上达到更高的要求。混合法是上述两种方法的结合。
（4）“实现问题”就是根据建立的数学模型和精度，采用某种数值计算方法，将模型方程转换为适合在计算机上运行的公式和方程，通过计算来使之正确的反映系统各变量动态性能，得到可靠的仿真结果。
（5）数值积分法应该遵循的原则是在满足系统精度的前提下，提高数值运算的速度和并保证计算结果的稳定。

2-2.用matlab语言求下列系统的状态方程、传递函数、零极点增益、和部分分式形式的模型参数，并分别写出其相应的数学模型表达式：
(1) G（s）=

(2) =
y=[0 2 0 2] X

（1） 解：（1）状态方程模型参数:编写matlab程序如下
>> num=[1 7 24 24];
>> den=[1 10 35 50 24];
>> [A B C D]=tf2ss(num,den)
得到结果：A=,B=,C=,D=[0]
所以模型为: =X+u,y=X
(2) 零极点增益:编写程序 >> num=[1 7 24 24];
>> den=[1 10 35 50 24];
>> [Z P K]=tf2zp(num,den)
得到结果Z= -2.7306 + 2.8531 , -2.7306 - 2.8531i ,-1.5388
P= -4, -3 ,-2 ,-1
K=1
(3) 部分分式形式:编写程序>> num=[1 7 24 24];
>> den=[1 10 35 50 24];
>> [R P H]=residue(num,den)
得到结果R= 4.0000 ,-6.0000, 2.0000, 1.0000
P= -4.0000, -3.0000 , -2.0000 ,-1.0000
H=[]
G(s)=

（2）解：（1）传递函数模型参数：编写程序>> A=[2.25 -5 -1.25 -0.5
2.25 -4.25 -1.25 -0.25
0.25 -0.5 -125 -1
1.25 -1.75 -0.25 -0.75]；
>> B=[4 2 2 0]';
>> C=[0 2 0 2];
>> D=[0];
>> [num den]=ss2tf(A,B,C,D)
得到结果
num = 0 4.0000 14.0000 22.0000 15.0000
den =1.0000 4.0000 6.2500 5.2500 2.2500

(2) 零极点增益模型参数:编写程序>> A=[2.25 -5 -1.25 -0.5
2.25 -4.25 -1.25 -0.25
0.25 -0.5 -1.25 -1
1.25 -1.75 -0.25 -0.75]；
>> B=[4 2 2 0]';
>> C=[0 2 0 2];
>> D=[0];
>> [Z,P,K]=ss2zp(A,B,C,D)

得到结果Z =-1.0000 + 1.2247i -1.0000 - 1.2247i -1.5000
P= -0.5000 + 0.8660i -0.5000 - 0.8660i -1.5000 -1.5000
K = 4.0000

表达式

（3）部分分式形式的模型参数：编写程序>> A[2.25 -5 -1.25 -0.5
2.25 -4.25 -1.25 -0.25
0.25 -0.5 -1.25 -1
1.25 -1.75 -0.25 -0.75]；
>> B=[4 2 2 0]';
>> C=[0 2 0 2];
>> D=[0];
>> [num den]=ss2tf(A,B,C,D)
>> [R,PH]=residue(num,den)
得到结果R = 4.0000 -0.0000 0.0000 - 2.3094i 0.0000 + 2.3094i
P = -1.5000 -1.5000 -0.5000 + 0.8660i -0.5000 - 0.8660i
H =[]

2-3.用欧拉法求下面系统的输出响应y(t)在0≤t≤1上，h=0.1时的数值。

要求保留4位小数，并将结果与真解比较。
解：欧拉法（前向欧拉法,可以自启动）其几何意义：把f(t,y)在[]区间内的曲边面积用矩形面积近似代替。利用matlab提供的m文件编程，得到算法公式。如下所示

（1） m文件程序为 h=0.1;
disp('函数的数值解为'); %显示 ‘’中间的文字%
disp('y='); %同上%
y=1;
for t=0:h:1
m=y;
disp(y); %显示y的当前值%
y=m-m*h;
end
保存文件q2.m
在matalb命令行中键入>> q2
得到结果 函数的数值解为
y= 1 0.9000 0.8100 0.7290 0.6561 0.5905 0.5314 0.4783 0.4305 0.3874 0.3487

（2）另建一个m 文件求解在t[0,1]的数值 （ %是的真解%）
程序为h=0.1;
disp('函数的离散时刻解为');
disp('y=');
for t=0:h:1
y=exp(-t);
disp(y);
end 保存文件q3.m
在matalb命令行中键入>> q3
函数的离散时刻解为
y= 1 0.9048 0.8187 0.7408 0.6703 0.6065 0.5488 0.4966 0.4493 0.4066 0.3679

比较欧拉方法求解与真值的差别

欧拉
1
0.9000
0.8100
0.7290
0.6561
0.5905
0.5314
0.4783
0.4305
0.3874
0.3487
真值
1
0.9048
0.8187
0.7408
0.6703
0.6065
0.5488
0.4966
0.4493
0.4066
0.3679
误差
0
-0.0048
-0.0007
–0.0118
–0.0142
–0.0160
–0.0174
–0.0183
–0.0188
-0.0192
-0.0192
显然误差与为同阶无穷小，欧拉法具有一阶计算精度，精度较低，但算法简单。

2-4用二阶龙格库塔法求解2-3的数值解，并于欧拉法求得的结果比较。
解：我们经常用到 预报-校正法 的二阶龙-格库塔法， 此方法可以自启动，具有二阶计算精度，几何意义：把f(t,y)在[]区间内的曲边面积用上下底为和、高为h的梯形面积近似代替。利用matlab提供的m文件编程，得到算法公式。如下所示
（1）m文件程序为 h=0.1;
disp('函数的数值解为');
disp('y=');
y=1;
for t0:h:1
disp(y);
k1=-y;
k2=-(y+k1*h);
y=y+(k1+k2)*h/2;
end
保存文件q4.m
在matlab的命令行中键入 >> q4 显示结果为
函数的数值解为
y= 1 0.9050 0.8190 0.7412 0.6708 0.6071 0.5494 0.4972 0.4500 0.4072 0.3685

2） 比较欧拉法与二阶龙格-库塔法求解.（误差为绝对值）

真值
1
0.9048
0.8187
0.7408
0.6703
0.6065
0.5488
0.4966
0.4493
0.4066
0.3679
龙库
1
0.9050
0.8190
0.7412
0.6708
0.6071
0.5494
0.4972
0.4500
0.4072
0.3685
误差
0
0.0002
0.0003
0.0004
0.0005
0.0006
0.0006
0.0006
0.0007
0.0006
0.0006

明显误差为得同阶无穷小，具有二阶计算精度，而欧拉法具有以阶计算精度，二阶龙格-库塔法比欧拉法计算精度高。

2-5．用四阶龙格-库塔法求解题2-3数值解，并与前两题结果相比较。
解：四阶龙格-库塔法表达式，其截断误差为同阶无穷小，当h步距取得较小时，误差是很小的.
(1)编辑m文件程序h=0.1;
disp('四阶龙格-库塔方法求解函数数值解为');
disp('y=');
y=1;
for t=0:h:1
disp(y);
k1=-y;
k2=-(y+k1*h/2);
k3=-(y+k2*h/2);
k4=-(y+k3*h);
y=y+(k1+2*k2+2*k3+k4)*h/6;
end 保存文件q5.m
在matlab命令行里键入>> q5
得到结果 四阶龙格-库塔方法求解函数数值解为
y= 1 0.9048 0.8187 0.7408 0.6703 0.6065 0.5488 0.4966 0.4493 0.4066 0.3679
(2)比较这几种方法：
对于四阶龙格-库塔方法

真值
1
0.9048
0.8187
0.7408
0.6703
0.6065
0.5488
0.4966
0.4493
0.4066
0.3679
龙库
1
0.9048
0.8187
0.7408
0.6703
0.6065
0.5488
0.4966
0.4493
0.4066
0.3679
误差
0
0
0
0
0
0
0
0
0
0
0

显然四阶龙格-库塔法求解精度很高，基本接近真值。三种方法比较可以得到
精度（四阶 ） 〉精度（二阶） 〉精度（欧拉）
2-6．已知二阶系统状态方程为
写出取计算步长为h时，该系统状态变量X=[]的四阶龙格-库塔法递推关系式。
解：四阶龙格-库塔法表达式
所以状态变量的递推公式可以写作：
A=,B=,可以写成
则递推形式

2-7单位反馈系统的开环传递函数已知如下

用matlab语句 、函数求取系统闭环零极点，并求取系统闭环状态方程的可控标准型实现。
解：已知开环传递函数，求得闭环传递函数为
在matlab命令行里键入>> a=[1 0];
>> b=[1 4.6];
>> c=[1 3.4 16.35];
>> d=conv(a,b)；
>> e=conv(d,c)
e = 1.0000 8.0000 31.9900 75.2100 0
>> f=[0 0 0 5 100];
>> g=e+f
g = 1.0000 8.0000 31.9900 80.2100 100.0000
%以上是计算闭环传递函数的特征多项式%
>> p=roots(g) %计算特征多项式的根，就是闭环传递函数的极点%
p =
-0.9987 + 3.0091i
-0.9987 - 3.0091i
-3.0013 + 0.9697i
-3.0013 - 0.9697i
>> m=[5 100];
>> z=roots(m)
z = -20 %计算零点%
综上：当闭环传函形如时，可控标准型为：
；
所以可控标准型是

2-8用matlab语言编制单变量系统三阶龙格-库塔法求解程序，程序入口要求能接收状态方程各系数阵（A,B,C,D）,和输入阶跃函数r(t)=R*1(t);程序出口应给出输出量y（t）的动态响应数值解序列
。
解：m文件为：function y=hs(A,B,C,D,R,T,h) %T为观测时间,h为计算步长，R为输入信号幅值%
disp('数值解为');
y=0;
r=R;
x=[0;0;0;0];
N=T/h;
for t=1:N;
k1=A*x+B*R;
k2=A*(x+h*k1/3)+BR;
k3=A*(x+2*h*k2/3)+B*R;
x=x+h*(k1+3*k3)/4;
y(t)=C*x+D*R;
end
在命令行里键入A= B= C= D= R= T= h=
y=hs(A,B,C,D,R,T,h) 得到结果。

2-9．用题2-8仿真程序求解题2-7系统的闭环输出响应y(t).
解：A=,B=,C=,D=[0]
在命令行里键入>> A=[0 1 0 0
0 0 1 0
0 0 0 1
-100 -80.21 -31.99 -8];
>> B=[0 0 0 1]';
>> C=[-100 5 0 0];
>> D=[0];
>> T=1;
>> R=1;
>> h=0.01;
>> y=hs(A,B,C,D,R,T,h)
数值解为
0

8.3333e-007

5.8659e-006

1.8115e-005

3.9384e-005

7.0346e-005
。
。
。
。 %仅取一部分%
2-10.用式（2-34）梯形法求解试验方程，分析对计算步长h有何限制，说明h对数值稳定性的影响。
解：编写梯形法程序为
得到 稳定系统最终渐进收敛。
系统稳定则 计算得。
h的选取不能超出上述范围，否则系统不稳定。

2-11如图2-27所示斜梁滚球系统，若要研究滚球在梁上的位置可控性，需首先建立其数学模型，已知力矩电机的输出转矩M与其电流i成正比，横梁为均匀可自平衡梁（即当电机不通电且无滚球时，横梁可处于=0的水平状态），是建立系统的数学模型，并给出简化后系统的动态结构图。
解：设球的质心到杆的距离为0，该系统为特殊情况下的球棒系统。另令分别表示棒的惯量、球的质量和球的惯量。则球质心的位置和速度为

其中，。因而动能的移动部分为
因而动能的移动部分为
球棒系统的旋转动能为
因而，系统总的动能等于

其中为常数。
此系统的拉格朗日方程组为

综合以上公式的系统的方程组为

设系统在平衡点附近，，，则系统方程可化为

对上式进行拉普拉斯变换并化简后可得到。
参考文献：
[1] Hauser, S. Sestry, and P. Kokotovic. “Nonlinear control via approximate input-output linearization”. IEEE Trans. on Automatic Control, vol.37:pp.392-398, 1992.
[2] R. Sepulchre. “Slow peaking and low-gain designs for global stabilization of nonlinear systems”. submitted for IEEE TAC 1999.
[3] R. Sepulchre, M. Jankovic, and P. Kokotovic Constructive Nonlinear Control. Springer-Verlag, 1997.
[4] R. Teel. “Using Saturation to stabilize a class of single-input partially linear composite systems”. IFAC NOLCOS'92 Symposium, pages 369-374, June 1992.

2-12如图2-28所示双水箱系统中，为流入水箱1的液体流量，为流出水箱2的液体流量，试依据液容与液阻的概念，建立的系统动态结构图。
解：根据液容和液阻的概念，可分别列出两个水箱的数学模型

对上式进行在零初始条件下进行拉普拉斯变换得

化简后可得

第三章 习题
4-2设典型闭环结构控制系统如图4-47所示，当阶跃输入幅值 时，用sp4_1.m求取输出的响应。

解：用sp4_1.m求解过程如下：
在MATLAB语言环境下，输入以下命令语句
>> a=[0.016 0.864 3.27 3.42 1];
>> b=[30 25];
>> X0=[0 0 0 0]; %系统状态向量初值为零
>> V=2; %反馈系数
>> n=4;
>> T0=0;Tf=10;
>> h=0.01;R=20 ; %仿真步长h=0.01，阶跃输入幅值
>> sp4_1 %调用sp4_1.m函数
>> plot(t,y)
运行结果为：

附：sp4_1.m函数为
b=b/a(1);a=a/a(1);A=a(2:n+1);
A=[rot90(rot90(eye(n-1,n)));-fliplr(A)];
B=[zeros(1,n-1),1]';
m1=length(b);
C=[fliplr(b),zeros(1,n-m1)];
Ab=A-B*C*V;
X=X0';
y=0;t=T0;
N=round((Tf-T0)/h);
for i=1:N
K1=Ab*X+B*R;
K2=Ab*(X+h*K1/2)+B*R;
K3=Ab*(X+h*K2/2)+B*R;
K4=AbX+h*K3)+B*R;
X=X+h*(K1+2*K2+2*K3+K4)/6;
y=[y,C*X];
t=[t,t(i)+h];
end
4-4系统结构图如图4-48，写出该系统的联结矩阵和，并写出联结矩阵非零元素阵。

解：根据图4-48中,拓扑连结关系，可写出每个环节输入受哪些环节输出的影响，
现列如入下:

把环节之间的关系和环节与参考输入的关系分别用矩阵表示出来，

即=，=，

4-6若系统为图4-5b 双输入-双输出结构，试写出该系统的联接矩阵，，说明应注意什么？

解：根据图4-5b中,拓扑连结关系，可列写如下关系式:

转换成矩阵形式为

所以联接矩阵=，=
此时应注意输入联接矩阵变为型。
4-8求图4-49非线性系统的输出响应y(t),并与无非线性环节情况进行比较。

解：（1）不考虑非线性环节影响时，求解过程如下：
1） 先将环节编号标入图中。
2) 在MATLAB命令窗口下，按编号依次将环节参数输入P阵；
>> P=[0.1 1 0.5 1;0 1 20 0;2 1 1 0;10 1 1 0];
3) 按各环节相对位置和联接关系，有联接矩阵如下：
， ，所以非零元素矩阵
>> WIJ=[1 0 1;1 4 -1;2 1 1;3 2 1;4 3 1] ;
4）由于不考虑非线性影响，则非线性标志向量和参数向量均应赋零值；
>> Z=[0 0 0 0];S=[0 0 0 0];
5）输入运行参数：开环截至频率约为1，故计算步长h取经验公式值，即
，取h=0.01；每0.25秒输出一点。故取=25。
>>h=0.01;
>>L1=25;
>>n=4;
>>T0=0
>>Tf=20;
>>nout=4;
>>Y0=10;
>>sp4_4;
>> plot(t,y,'r')
>> hold on
运行结果如图中红色实线所示。
(2)考虑非线性环节N影响时，只需将非线性标志向量Z和参数向量S的相应分量正确输入即可。
在MATLAB命令窗口中输入下列语句：
>> Z=[4 0 0 0];S=[5 0 0 0]; %第一个线性环节后有饱和非线性，参数值为5。
>> sp4_4;
>> plot(t,y,'--')
运行结果如图中蓝色虚线所示。

从图中可以清楚的地看出，饱和非线性环节对线性系统输出响应的影响。
附：sp4_4函数为:
A=P(:,1);B=P(:,2);
C=P(:,3);D=P(:,4);
m=length(WIJ(:,1));
W0=zeros(n,1);W=zeros(n,n);
for k=1:m
if (WIJ(k,2)==0); W0(WIJ(k,1))=WIJ(k,3);
else W(WIJ(k,1),WIJ(k,2))=WIJ(k,3);
end;
end;
for i=1:n
if(A(i)==0);
FI(i)=1;
FIM(i)=h*C(i)/B(i);
FIJ(i)=h*h*(C(i)/B(i))/2;
FIC(i)=1;FID(i)=0;
if(D(i)~=0);
FID(i)=D(i)/B(i);
else
end
else
FI(i)=exp(-h*A(i)/B(i));
FIM(i)=(1-FI(i))*C(i)/A(i);
FIJ(i)=h*C(i)/A(i)-FIM(i)*B(i)/A(i);
FIC(i)=1;FID(i)=0;
if(D(i)~=0);
FIC(i)=C(i)/D(i)-A(i)/B(i);
FID(i)=D(i)/B(i);
else
end
end
end
Y=zeros(n,1);X=Y;y=0;
Uk=zeros(n,1);Ubb=Uk;
t=T0:h*L1:Tf;N=length(t);
for k=1:N-1
for i=1:L1
Ub=Uk;
Uk=W*Y+W0*Y0;
for i=1:n
if(Z(i)~=0)
if (Z(i)==1)
Uk(i)=satu(Uk(i),S(i));
end
if(Z(i)==2)
Uk(i)=dead(Uk(i),S(i));
end
if(Z(i)==3)
[Uk(i),Ubb(i)]=backlash(Ubb(i),Uk(i),Ub(i),S(i));
end
end
end
Udot=(Uk-Ub)/h;
Uf=2*Uk-Ub;
X=FI'.*X+FIM'.*Uk+FIJ'.*Udot;
Yb=Y;
Y=FIC'.*X+FID'.*Uf;
for i=1:n
if(Z(i)~=0)
if (Z(i)==4)
Y(i)=satu(Y(i),S(i));
end
if(Z(i)==5)
Y(i)=dead(Y(i),S(i));
end
if(Z(i)==6)
[Y(i),Ubb(i)]=backlash(Ubb(i),Y(i),Yb(i),S(i));
end
end
end
end
y=[y,Y(nout)];
end
附：饱和非线性函数satu.m 为：
function Uc=satu(Ur,S1)
if(abs(Ur)>=S1)
if(Ur>0)
Uc=S1;
else Uc=-S1;
end
else Uc=Ur;
end

4-10采样控制系统如图4-50所示，编写程序实现对该系统的仿真分析。（提示：连续部分按环节离散化考虑）图中，为典型数字PID控制器；=0.65为比例系数；=0.7为积分时间常数；=0.2为微分时间常数；
为具有纯滞后特性的典型二阶控制对象；；；。
解：在控制对象前引入零阶保持器，将连续环节部分按环节离散化：
==
设，为简化运算及编程，取为T的整数倍

=
对上式进行Z逆变换，得到

由此可编写仿真程序。
在MATLAB命令窗口中输入下列语句：
>> KP=0.65;TI=0.7;TD=0.2;
>> T1=0.3;a=1/T1;T3=0.4;
>> T=0.1;h=0.001;Tf=10;
>>hh

编写M脚本文件，存为hh.m。
%离散化后各参数为：
A=1-a*h*exp(-ah)-exp(-a*h);B=exp(-2*a*h)-exp(-a*h)+a*h*exp(-a*h);C=2*exp(-a*h);D=exp(-2*a*h);
P=KP*(1+T/TI+TD/T);H=KP*(1+2*TD/T);M=KP*TD/T;
%系统初始值为：
E=zeros(1,3);U=zeros(1,2+T3/T+1);Y=zeros(1,2+T3/h+1);R=1;
yk=0;yt=0;t=0;
%仿真迭代运算：
for K1=1:Tf/T
ek=R-Y(1);
E=[ek,E(1:2)];
uk=P*E(1)-H*E(2)+M*E(3)+U(1);
U=[uk,U(1:(2+T3/T))];
for K2=1:T/h
yk=A*U(T3/T+1+1)+B*U(T3/T+2+1)+C*Y(T3/h+1)-D*Y(T3/h+2);
Y=[yk,Y(1:(2+T3/h))];
end
yt=[yt,yk];
t=[t,K1*T];
end
%输出波形：
plot(t,yt)

运行结果为：

此题可以用SIMULINK仿真进行验证：
建立SIMULINK仿真模型：

运行结果为：

第四章习题
5-1 设控制系统的开环传递函数为

试画出该系统的根轨迹。

解: 在Matlab窗口中输入下列命令：
num=[1 1];
a=[1 0];
b=[1 -1];
c=[1 4 16];
d=conv(a,b);
den=conv(d,c);
rlocus(num,den)
grid on
可得到系统的根轨迹如下图所示：

5-2 某反馈控制系统的开环传递函数为

试绘制其根轨迹。

解：在MATLAB命令窗口中输入下列命令：
num=1;
den=conv(conv([1,0],[1,4]),[1,4,2