《控制系统数字仿真与CAD》张晓华版课后答案doc

2-2.用MATLAB语言求下列系统的状态方程
（1）状态方程模型参数:编写matlab 程序如下
>>num=[1 7 24 24];
>>den=[1 10 35 50 24];
>>[AB C D]=tf2ss(num,den)
(2) 零极点增益:编写程序>> num=[1 7 24 24];
>> den=[1 10 35 50 24];
>> [ZPK]=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 PH]=residue(num,den)
得到结果R= 4.0000 ,-6.0000, 2.0000, 1.0000
P= -4.0000, -3.0000 , -2.0000 ,-1.0000
H=[]
（2）解：（1）传递函数模型参数：
编写程序>>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)
得到结果
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)
（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,P,H]=residue(num,den)
2-6
(1)解：m 文件程序为h=0.1;
disp('函数的数值解为'); %显示‘’中间的文字%
disp('y='); %同上%
y=1;
fort=0:h:1
m=y;
disp(y); %显示y 的当前值%
y=m-m*h;
end
保存文件da2.m
在matalb 命令行中键入>> da
得到结果函数的数值解为
y= 1 0.9000 0.8100 0.7290 0.6561 0.5905 0.5314 0.4783 0.4305 0.3874 0.3487 （2）另建一个m文件求解y=e在t∈[0,1]的数值
程序为h=0.1;
disp('函数的离散时刻解为');
disp('y=');
fort=0:h:1
y=exp(-t);
disp(y);
end 保存文件da3.m
在matalb 命令行中键入>> da3
函数的离散时刻解为
y= 1 0.9048 0.8187 0.7408 0.6703 0.6065 0.5488 0.4966 0.4493 0.4066 0.3679 2-5.
解：
(1)编辑m 文件程序h=0.1;
disp('四阶龙格-库塔方法求解函数数值解为');
disp('y=');
y=1;
fort=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-7
>>a=[1 0];
>> b=[14.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=[00 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 %计算零点计算零点计算零点计算零点%
3-2．进行下列计算，给出不使用for 和while 等循环语句的计算方法。

（1）6302
i i k ==∑
解:根据等比数列求和方法，再利用matlab 中的m 文件，编写程序求解。

M 文件为 n=64;
q=2;
k=(1-q^n)/(1-q);
disp('k 的值为');
disp(k);
保存文件q1.m
在matlab 命令框中输入 >> q1
k 的值为
1.8447e+019
(2)求出y=x*sin(x) 在0<x<100条件下的每个峰值
解：画出图形
>> x=0:0.01:100;
>> y=x.*sin(x);
>> plot(x,y);
>> grid on
>> title('y=x*sin(x)')
>> xlabel('x')
>>ylabel('y')
从图形中不难看出峰值点取决于函数sin(x),即在sin(x)为峰值时，y 就得到峰值。

所以求取函数的峰值转化为求取正弦函数波峰问题。

而sin(x)在x=
2
π+2k π(k 为整数)，所以求取y 在上述x 时刻的数值就是峰值。

在matlab 命令行里键入
>> x=pi/2:pi*2:100；
>> y=x.*sin(x)
得到结果
y=1.5708 7.8598 14.1481 20.4350 26.7198 33.0019 39.2804 45.5549 51.8245
58.0887 64.3467 70.5978 76.8414 83.0769 89.303 95.5204
3-8
解：：：（1）不考虑非线性环节影响时，求解过程如下：
1）先将环节编号标入图中。

2) 在 MATLAB 命令窗口下，按编号依次将环节参数输入 P 阵；
>> P=[0.1 10.5 1;01 200;2 11 0;1011 0];
3) 按各环节相对位置和联接关系，有联接矩阵如下：
>>WIJ=[1 0 1;1 4 -1;2 1 1;3 2 1;4 3 1];
4）由于不考虑非线性影响，则非线性标志向量和参数向量均应赋零值；
>>Z=[0 0 00];S=[0 0 0 0];
5）输入运行参数：开环截至频率1C L ω约为1，故计算步长 h 取经验公式值，即
150c
h ω≤=0.02, 取 h=0.01；每 0.25秒输出一点。

故取L1 =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 的相应分量
正确输入即可。

在MA TLAB 命令窗口中输入下列语句：
>> Z=[4 00 0];S=[5 0 0 0]; %第一个线性环节后有饱和非线性，参数值为5。

>> sp4_4;
>> plot(t,y,'--')
运行结果如图中虚线所示。

附1：
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);
elseW(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
fori=1:L1
Ub=Uk;
Uk=W*Y+W0*Y0;
fori=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;
fori=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
附2：饱和非线性函数satu.m 为：functionUc=satu(Ur,S1)
if(abs(Ur)>=S1)
if(Ur>0)
Uc=S1; else Uc=-S1; end
else Uc=Ur; end。