Matlab作业答案.docx

控制系统的MATLAB仿真与设计 - 实验作业
第二章 数据结构及其运算
2.1实验程序：
x=[12 23 34 45 56 67 78] %逐个输入数组中的元素值
x(6) %读数组中的第6个元素
x([1 3 5]) %返回数组的第1、3、5个元素构成的子数组
x(4:end) %返回第4个元素与其后所有元素构成的子数组
x(find(x>70)) %返回数组中大于70的元素的值
MATLAB分步截屏：
创建一个有7个元素的一维数组：

直接寻访一维数组的第6个元素：

寻访一维数组的第1，3，5个元素：

(3)寻访一维数组中第4个至最后一个元素：

寻访一维数组中大于70 的元素：

2.3 实验程序：
A=zeros(2,5) %生成2行5列的全0数组
A(1:10)=-4:5 %将单下标为1到10的元素分别赋值为-4到5
L=abs(A)>3 %判断其绝对值是否大于3，大于则L对应的位置为值1，否则为L值为0
islogical(L) %判断L是否为逻辑数组，是逻辑数组则返回1，否则返回0
X=A(L) % X为逻辑数组L对应为1的位置的数据的列向量
MATLAB分步截屏：
生成数组：

寻访数组A中绝对值大于3的元素：

输出数组A中绝对值大于3的元素：

2.5实验程序：
p1=[1,0,2]; %定义多项式p1(s)=s2+1
p2=[1,4]; %定义多项式p2(s)=s+4
p3=[1,1]; %定义多项式p3(s)=s+1
p4=conv(p1,conv(p2,p3)); %定义多项式p4(s)=p1(s)p2(s)p3(s)
p5=[1 0 1 1]; %定义多项式p5(s)=s3+s+1
[q,r]=deconv(p4,p5); % 求p4(s)/p5(s)
disp(['商:',poly2str(q,'s')]);
disp(['余式:',poly2str(r,'s')])
MATLAB实验截屏：

2.6实验程序：
A=[12 23 34;45 56 67;78 89 90]; %定义3阶数组A
p1A=poly(A); %创建方阵A的特征多项式p1A
p2A=poly2str(p1A,'s') %将拟合后的p1A多项式系数转换为字符形式的特征多项式函数p2A
MATLAB实验截屏：

第四章 M文件和程序设计
4.2实验程序exp4_2：
x=input('请输入半径r的值：'); %将用户键入量赋给x
s=f1(x); %调用函数f1(x)计算面积S
c=f2(x); %调用函数f2(x)计算面积C
disp(['圆的面积S为:',num2str(s)]) %用num2str将数值变量转换为字符变量并用disp输出带数值的文本信息
disp(['圆的周长C为:',num2str(c)])
实验程序-面积函数文件f1：
function s=f1(x) %自定义函数f1(x)
s=pi*x*x;
end
实验程序-周长函数文件f2：
function c=f2(x) %自定义函数f2(x)
c=2*pi*x;
end
MATLAB截屏：

4.3实验程序exp4_3：
sum=0; %定义变量sum
n=input('请输入n的值：') %并将用户键入量赋给变量n，并显示“请输入n的值”
for i=1:n % i从1取到n
sum=sum+1/i^2; %依次求sum的值
end
disp(['计算结果为：',num2str(sum)]) %用num2str将数值变量转换为字符变量并用disp输出带数值的文本信息

MATLAB截屏：

4.5实验程序exp4_5：
sum=0; %变量sum的初值是0
n=0; %变量n的初值是0
while sum<2000 %当sum值小于2000时
n=n+1; %n=n+1
sum=sum+n; %sum为n的累加和
end
disp(['最大的n值为：',num2str(n-1)]) %最后一次循环，sum+n后得到的sum值大于2000，因此最后输出的最大n值应为n-1
MATLAB截屏：

4.6(1)使用for循环时：
实验程序exp4_6_1：
sum=0; %定义变量sum，其初值为0
for i=0:63 % i从0取到63
sum=sum+2^i; %依次求sum的累加
end
disp(['K=',num2str(sum)])
MATLAB截屏：

(2)使用while循环时：
实验程序exp4_6_2：
sum=0; %定义变量sum，其初值为0
i=0; %定义变量i，其初值为0
while i<=63 %当i≦63时
sum=sum+2^i; %依次求sum的累加
i=i+1;
end
disp(['K=',num2str(sum)])
MATLAB截屏：

第五章 数值计算
5.1用LU分解求解下列线性方程组：
2x1+x2-5x3+x4=13x1-5x2+7x4=-92x1+x3-x4=6x1+6x2-x3-4x4=0
实验程序exp5_1：
A=[2,1,-5,1;1,-5,0,7;2,0,1,-1;1,6,-1,-4];
b=[13,-9,6,0]';
[L,U]=lu(A); %进行LU分解
x=U\(L\b) %线性方程组的解
MATLAB结果：

5.2设fx,y,z=x+y24x+z2y+2z ，求函数f在（0.5,0.5,0.5）附近的最小值：
实验程序exp5_2：
fun=inline('x(1)+x(2)^2/(4*x(1))+x(3)^2/x(2)+2/x(3)','x')
%使用内联函数构造函数fun
[x,fval]=fminsearch(fun,[0.5,0.5,0.5]) %求函数fun在（0.5,0.5,0.5）附近的最小值
MATLAB结果：

5.5有初值问题：
y'=y2-t-24(t+1) 0≤t≤1y0=2
试求其数值解，并与下式精确解作比较：yt=t+1+1
实验程序exp5_5：
OdeFun=inline('(y^2-t-2)/(4*(t+1))') %内建函数编写待求微分方程
ts=[0,1]; %自变量求解区间
y0=2; %初值条件
[t,y]=ode23(OdeFun,ts,y0); %2、3阶龙格库塔法求解微分方程数值解
y1=sqrt(t+1)+1; %微分方程精确解
[t,y,y1] %以三列矩阵显示，第一列为采样点t对应的点，第二列为期数值解，第三列为其精确解
MATLAB结果：

第七章 Simulink基础
7.1单位负反馈系统的开环传递函数为
Gs=1000s(0.1s+1)(0.001s+1)
应用Simulink仿真系统构建其阶跃响应曲线。
（1）Simulink仿真：

（2）系统输出：

（3）验证：
应用Matlab函数step进行验证如下，系统单位阶跃响应的稳态值为1：

7.2系统微分方程为x1'=x2tx2'=x2e-0.5t ，试建立系统Simulink仿真模型。
Simulink仿真：
①
②
③
④
⑤
⑥
⑦
⑧
⑨

仿真图中位置处的信号说明如下：
① 时钟信号t
② 信号x2
③ x1'=x2t
④ 输出信号x1
⑤ -0.5t
⑥ e-0.5t
⑦ 信号x2
⑧ x2'=x2e-0.5t
⑨ 信号x2
7.3利用Simulink构建函数曲线y=5t2+16。
（1）Simulink搭建仿真电路：
将时钟信号做平方数学运算得到单位抛物线函数，扩大5倍后与常数输入16叠加作为系统给定输入。

（2）系统输出：
①示波器Scope输出函数曲线 y=5t2+16：

②示波器Scope1输出单位抛物线曲线 y=t2：

（3）验证：
在命令窗口输入程序：>> x=0:0.001:10;
>> y=5.*x.*x+16;
>> plot(x,y,'k-')
>> grid on

与示波器Scope得出的仿真结果相同。
7.4利用Simulink仿真I=01xln(1+x)dx。
（1）Simulink截图：
将时钟信号与常数输入1相加后做自然对数ln数学函数运算，再与时钟信号相乘叠加后做积分运算，将仿真时间设为1s，输出。

（2）系统输出：
①数字显示模块显示积分结果为0.25：

②示波器Scope输出积分结果随积分时间t的变化：

③示波器Scope1输出积分函数曲线I=01xln(1+x)dx：

（3）验证：
①积分结果验证：
输入程序：>> x=0:0.001:1;
>> y=x.*log(1+x);
>> trapz(x,y)
s =

0.2500
积分结果为0.25，Simulink仿真在t=1s时输出为0.25，两者结果相同。
②积分函数曲线验证：
输入程序：>> plot(x,y)
>> plot(x,y,'k-')
>> grid on
输出积分函数图像如下，与Simulink仿真示波器Scope1输出结果相同：

7.6利用使能子系统构成一个正弦半波整流器。
（1）Simulink仿真：
使用电平触发的Enabled Subsystem子系统，当使能信号为高电平的时候运行子系统，得出正弦半波整流波形。

（2）系统输出：

第八章 控制系统的数学模型
8.1创建连续二阶系统的传递函数模型
(1)G(s)=5s2+2s+2
实验程序exp8_1_1：
num=[0,0,5]; %分子多项式的系数
den=[1,2,2]; %分母多项式的系数
sys=tf(num,den) %建立系统的TF模型
MATLAB截屏：

(2)G(s)=5s2+2s+2e-2s
实验程序exp8_1_2：
num=[0,0,5];
den=[1,2,2];
sys=tf(num,den);
sys.inputdelay=2 %输入信号的滞后时间为2秒
MATLAB截屏：

8.2已知系统的传递函数为G(s)=2(s+0.5)(s+0.1)2+1，建立系统的传递函数模型，并转换为零极点模型
实验程序exp8_2：
s=tf('s');
G1=2*(s+0.5)/((s+0.1)^2+1) %直接用数学运算符建立系统的TF模型并显示
num=G1.num{1}; %通过属性名访问分子的系数并赋给数组num
den=G1.den{1}; %通过属性名访问分母的系数并赋给数组den
sys=tf(num,den); %以num和den建立TF模型
[z,p,k]=tf2zp(num,den); %将TF对象的属性转换为ZPK对象的属性
G2=zpk(z,p,k) %建立系统的ZPK模型
MATLAB截屏：

8.4已知系统的方框图如下，其中R1=1，R2=2，C1=3，C2=4，计算系统的∅s=C(s)R(s)。
实验程序exp8_4：
s=tf('s'); %直接用数学运算符建立各个环节的TF模型
sys1=1+0*s;
sys2=1/(3*s);
sys3=1/2+0*s;
sys4=1/(4*s);
sys5=-1+0*s;
sys6=-1+0*s;
sys7=-1+0*s;
sys8=1+0*s; %sys8是sys2、sys3之间银触点和汇合点之间的连线，因为前后信号不同，所以不能合并，需要建立成一个环节
sys=append(sys1,sys2,sys3,sys4,sys5,sys6,sys7,sys8)；%扩展成无连接状态空间模型
Q=[1 5 0;2 1 7;3 8 6; 4 3 0;5 2 0;6 4 0;7 3 0;8 2 0]; %写出系统的连接矩阵
inputs=1; %sys1的输入为系统的总输入
outputs=4; %sys4的输出为系统的总输出
sysc=connect(sys,Q,inputs,outputs) %生成组合后的整体系统

MATLAB截屏：

第九章 控制系统的分析
9.1实验程序：
num=[2,5,1];
den=[1,2,3];
sys=tf(num,den); %建立系统的传递函数模型
rlocus(sys); %绘制根轨迹图
figure %绘图命令
bode(num,den); %绘制伯德图
grid on
figure
nyquist(num,den); %绘制奈奎斯特图
MATLAB结果：
(1)根轨迹图：

(2)伯德图：

(3)奈奎斯特图：

9.2 实验程序：
num=[5,25,30];
den=[1,6,10,8];
sys=tf(num,den); %建立系统的传递函数模型
step(sys) %绘制系统的单位阶跃响应曲线
figure
impulse(sys) %绘制系统的理想单位脉冲响应曲线
MATLAB结果：
(1)单位阶跃响应：

(2)单位脉冲响应：

第十章 基于传递函数模型的控制系统设计
10.4考虑一个单位负反馈控制系统，其前向通道传递函数为
G0s=1s2(s+5)
试应用伯德图法设计一个超前校正装置Gcs=Kcα(Ts+1αTs+1) ，使得校正后系统的相角裕量γ=50°，幅值裕度Kg≥10dB，带宽ωb=1~2rad/s。其中，0<α<1。试问已校正系统的谐振峰值Mr和谐振角频率ωr的值各为多少？
（1）源程序exp10_4：
ng0=[1];
dg0=[1,5,0,0];
t=[0:0.01:5];
w=logspace(-3,2);
disp('校正前系统开环传递函数为：')
g0=tf(ng0,dg0)
b1=feedback(g0,1); %校正前系统闭环传函
[gm,pm,wcg,wcp]=margin(g0); %校正前系统频域参数
disp(['校正前系统频域参数:','相角裕度',num2str(pm),' ','幅值裕度',num2str(gm)])
Pm=55; %校正后系统的相角裕度（保留一定预留量）
[ng1,dg1]=fg_lead_pm(ng0,dg0,Pm,w); %利用超前校正进行校正
disp('校正环节传递函数为：')
g1=tf(ng1,dg1) %校正环节传递函数
disp('校正后系统开环传递函数为：')
g2=g0*g1 %校正后前向通道传函
[gm1,pm1,wcg1,wcp1]=margin(g2); %校正后系统参数
disp(['校正后系统频域参数:','相角裕度',num2str(pm1),' ','幅值裕度',num2str(gm1)])
bode(g0,'r--',g1,'b--',g2,'g',w)
grid on
b2=feedback(g2,1); %校正后系统闭环传函
[mag,phase,w]=bode(b2);%对校正后系统闭环传函bode图进行离散化
Mr=max(mag); %求取Mr
disp(['校正后系统谐振峰值Mr:',num2str(Mr)])
a=find(mag==Mr);%求Mr对应的脚标
wr=w(a); %求wr
disp(['校正后系统谐振峰值角频率wr:',num2str(wr)])
b=find(mag<=0.707*mag(1));%求幅值小于0.707倍零频幅值的脚标所组成的数组
wb=w(b(1));%求wb
disp(['校正后系统的带宽wb:',num2str(wb)])
源程序fg_lead_pm：
function [ngc,dgc]=fg_lead_pm(ng0,dg0,Pm,w)
[mu,pu]=bode(ng0,dg0,w);
[gm,pm,wcg,wcp]=margin(mu,pu,w);
alf=ceil(Pm-pm+5);
phi=(alf)*pi/180;
a=(1+sin(phi))/(1-sin(phi));
a1=1/a;
dbmu=20*log10(mu);
mm=-10*log10(a);
wgc=spline(dbmu,w,mm);
T=1/(wgc*sqrt(a));
ngc=[a*T,1];
dgc=[T,1];

（2）MATLAB仿真结果：

从仿真结果可知，校正后系统满足预期的相角欲量和带宽要求，但幅值裕量为8.8907，略微低于题目要求。此时通过增加相角预留量重新进行系统校正发现，校正后系统幅值裕度增加缓慢，但相较于量和带宽有明显提高（如下表所示）。为了保证系统的平稳性，选择预留量为ε=5°，此时相角裕度如仿真结果。

相角裕度
幅值裕度
带宽
相角预留量ε=0°
51.098
8.3933
1.7055
相角预留量ε=5°
55.1177
8.8907
1.9373
相角预留量ε=8°
57.3808
9.3199
1.9816
表10-5 选择不同的ε时系统的频域参数

10.5 考虑一个单位负反馈控制系统，其前向通道传递函数为
G0s=Ks(s+1)(s+4)
试应用伯德图法设计一个校正装置Gcs，使得校正后系统的静态速度误差常数Kv=10s-1，相角裕量γ=50°，幅值裕度Kg≥10dB。
（1）源程序exp10_5：
kk=40; %校正后系统开环增益kk=4*kv
ng0=kk*[1];
dg0=conv([1,0],conv([1,1],[1,4]));
t=[0:0.01:5];
w=logspace(-3,2);
disp('校正前系统开环传递函数为：')
g0=tf(ng0,dg0) %校正前系统开关传递函数
[gm,pm,wcg,wcp]=margin(g0); %校正前系统开环传递函数的频域特性参数
disp(['校正前系统频域参数:','相角裕度',num2str(pm),' ','幅值裕度',num2str(gm)])
Pm=64.5; %期望相角裕量（留有一定的预留量）
[ng1,dg1]=fg_lead_pm(ng0,dg0,Pm,w); %利用超前校正进行校正
disp('校正环节传递函数为：')
g1=tf(ng1,dg1) %校正环节传递函数
disp('校正后系统开环传递函数为：')
g2=g1*g0 %校正后前向通道传函
[gm1,pm1,wcg1,wcp1]=margin(g2); %校正后系统频域参数
disp(['校正后系统频域参数:','相角裕度',num2str(pm1),' ','幅值裕度',num2str(gm1)])
bode(g0,'b',g1,'r',g2,'g',w)
grid on
（2）MATLAB仿真结果：

第十一章 基于状态空间模型的控制系统设计
11.1考虑单输入系统：x=Ax+bu，其中：A=010001-1-5-6，b=001。利用状态反馈控制u=-Kx，希望该系统的闭环极点为s1,2=-2±j4和s3=-10。确定状态反馈增益矩阵K。
（1）源程序exp11_1：
%利用Bass-Gura算法配置单输入系统极点
A=[0 1 0;0 0 1;-1 -5 -6];
b=[0;0;1];
p=[-2+i*4;-2-i*4;-10];
disp('状态反馈增益矩阵K为：')
K=bass_pp(A,b,p)
（2）MATLAB仿真结果：

11.8 考虑系统：x=Ax+buy=cx，其中：A=020.610，b=01，c=01。试设计一个全维状态观测器，使得观测器所期望的特征值为：s1,2=-1.8±j2.4。
（1）源程序exp11_8：
A=[0 20.06;1 0];
b=[0;1];
c=[0 1];
D=0;
p=eig(A)'
p1=p;
p1(1:2)=[-1.8+i*2.4;-1.8-i*2.4];
%观测器极点配置
L=place(A',c',p1)'
eig(A-L*c)'
%观测器的状态观测值输出
[xh,x,t]=simobsv(A,b,c,D,L);

（2）MATLAB仿真结果：

26