Chapter11.3：MATLAB_SIMULINK在离散系统中的应用

该系列博客主要讲述Matlab软件在自动控制方面的应用，如无自动控制理论基础，请先学习自动控制系列博文，该系列博客不再详细讲解自动控制理论知识。
自动控制理论基础相关链接：https://blog.csdn.net/qq_39032096/category_10287468.html?spm=1001.2014.3001.5482
博客参考书籍：《MATLAB/Simulink与控制系统仿真》。

---

---

3.MATLAB/SIMULINK在离散系统中的应用

3.1 MATLAB/SIMULINK在离散系统中常用函数

【连续系统模型与离散系统模型转换函数】

函数	调用格式	函数说明
$c2d$	$\mathrm{sysd}=c2d(\mathrm{sysc},\mathrm{Ts}{,}^{\prime }\mathrm{metho}{\mathrm{d}}^{\prime })$	连续时间$\mathrm{LTI}$系统模型转换为离散时间系统模型
$d2c$	$\mathrm{sysc}=d2c(\mathrm{sysd}{,}^{\prime }\mathrm{metho}{\mathrm{d}}^{\prime })$	离散时间$\mathrm{LTI}$系统模型转换为连续时间系统模型
$d2d$	$\mathrm{sys}=d2d(\mathrm{sysd},\mathrm{Ts})$	离散时间系统模型转换为新的$\mathrm{Ts}$离散时间系统模型

【$\mathrm{method}$功能说明】

选项	功能说明	选项	功能说明
${}^{\prime }\mathrm{zo}{\mathrm{h}}^{\prime }$	对输入信号加零阶保持器	${}^{\prime }\mathrm{tusti}{\mathrm{n}}^{\prime }$	双线性变换方法
${}^{\prime }\mathrm{fo}{\mathrm{h}}^{\prime }$	对输入信号加一阶保持器	${}^{\prime }\mathrm{prewar}{\mathrm{p}}^{\prime }$	预先转折变换方法，即改进的双线性变换方法
${}^{\prime }\mathrm{im}{\mathrm{p}}^{\prime }$	脉冲不变变换方法	${}^{\prime }\mathrm{matche}{\mathrm{d}}^{\prime }$	零极点匹配变换方法

【离散系统时域响应函数】

函数名	调用格式	功能说明
$\mathrm{dstep}$	$\begin{array}{rl}& \mathrm{dstep}(\mathrm{dnum},\mathrm{dden},\mathrm{n})\\ & \mathrm{y}=\mathrm{dstep}(\mathrm{dnum},\mathrm{dden},\mathrm{n})\end{array}$	求离散系统单位阶跃响应
$\mathrm{dimpulse}$	$\begin{array}{rl}& \mathrm{dimpulse}(\mathrm{dnum},\mathrm{dden},\mathrm{n})\\ & \mathrm{y}=\mathrm{dimpulse}(\mathrm{dnum},\mathrm{dden},\mathrm{n})\end{array}$	求离散系统单位脉冲响应
$\mathrm{dlsim}$	$\begin{array}{rl}& \mathrm{dlsim}(\mathrm{dnum},\mathrm{dden},\mathrm{u})\\ & \mathrm{y}=\mathrm{dlsim}(\mathrm{dnum},\mathrm{dden},\mathrm{u})\end{array}$	求离散系统在输入$\mathrm{u}$下的响应

注：$\mathrm{n}$为采样次数，$\mathrm{u}$为输入函数；

【离散系统频域响应函数】

函数名	调用格式	功能说明
$\mathrm{dbode}$	$\begin{array}{rl}& \mathrm{dbode}(\mathrm{dnum},\mathrm{dden},\mathrm{Ts},\mathrm{w})\\ & [\mathrm{mag},\mathrm{phase},\mathrm{w}]=\mathrm{dbode}(\mathrm{dnum},\mathrm{dden},\mathrm{Ts},\mathrm{w})\end{array}$	离散$\mathrm{Bode}$图
$\mathrm{dnyquist}$	$\begin{array}{rl}& \mathrm{dnyquist}(\mathrm{dnum},\mathrm{dden},\mathrm{Ts},\mathrm{w})\\ & [\mathrm{re},\mathrm{im},\mathrm{w}]=\mathrm{dnyquist}(\mathrm{dnum},\mathrm{dden},\mathrm{Ts},\mathrm{w})\end{array}$	离散$\mathrm{Nyquist}$图
$\mathrm{dnichols}$	$\begin{array}{rl}& \mathrm{dnichols}(\mathrm{dnum},\mathrm{dden},\mathrm{Ts},\mathrm{w})\\ & [\mathrm{re},\mathrm{im},\mathrm{w}]=\mathrm{dnichols}(\mathrm{dnum},\mathrm{dden},\mathrm{Ts},\mathrm{w})\end{array}$	离散$\mathrm{Nichols}$图
$\mathrm{margin}$	$\begin{array}{rl}& \mathrm{margin}(\mathrm{dsys})\\ & [\mathrm{Gm},\mathrm{Pm},\mathrm{Wcg},\mathrm{Wcp}]=\mathrm{margin}(\mathrm{dsys})\end{array}$	离散$\mathrm{Bode}$图，显示频域性能参数

注：$\mathrm{Ts}$为采样周期，$\mathrm{mag}$为幅值向量，$\mathrm{phase}$为相交向量，$\mathrm{Gm}$为增益裕量，$\mathrm{Pm}$为相角裕量，$\mathrm{re}$为$\mathrm{Nyquist}$图或$\mathrm{Nichols}$图实部向量，$\mathrm{im}$为$\mathrm{Nyquist}$图或$\mathrm{Nichols}$图虚部向量；

3.2 实战部分

3.2.1 实战1

实验要求：已知一个连续线性系统如下图所示，其中：$G_p(s)=\frac{1}{s(s+1)}$，用零阶保持器方法、一阶保持器方法、双线性变换方法和根匹配方法将此连续系统离散化，其中采样周期为：$T_s=0.1\mathrm{s}$。

解：

% 实例Chapter11.3.2.1
clc;clear;num=[1];den=[1,1,0];G=tf(num,den);      % 连续系统传递函数模型
Ts=0.1;Gd1=c2d(G,Ts,'zoh');                    % 零阶保持器方法
Gd2=c2d(G,Ts,'foh');                    % 一阶保持器方法
Gd3=c2d(G,Ts,'tustin');                 % 双线性变换方法
Gd4=c2d(G,Ts,'matched');                % 零极点匹配方法Gd1,Gd2,Gd3,Gd4

% 结果显示：% 零阶保持器方法
Gd1 =0.004837 z + 0.004679----------------------z^2 - 1.905 z + 0.9048% 一阶保持器方法
Gd2 =0.001626 z^2 + 0.006344 z + 0.001547------------------------------------z^2 - 1.905 z + 0.9048% 双线性变换方法
Gd3 =0.002381 z^2 + 0.004762 z + 0.002381------------------------------------z^2 - 1.905 z + 0.9048% 零极点匹配变换方法
Gd4 =0.004761 z + 0.004761----------------------z^2 - 1.905 z + 0.9048

3.2.2 实战2

实验要求：已知一个连续系统如下图所示，其中：$G_1(s)=\frac{2}{s(s+30)}，G_2(s)=\frac{10}{s^2+6s+5}$，采样周期$T_s=0.1\mathrm{s}$，求系统的脉冲闭环传递函数。

解：

% 实例Chapter11.3.2.2
clc;clear;Ts=0.1;
num1=[2];den1=[1,30,0];
num2=[10];den2=[1,6,5];G1=tf(num1,den1);
G2=tf(num2,den2);% 采用零阶保持器方法进行系统变换
G1d=c2d(G1,Ts);
G2d=c2d(G2,Ts);Gd=G1d*G2d;
GHd=feedback(Gd,1);      % 闭环系统模型% 模型显示
G1d,G2d,GHd

% 结果显示：% G1的离散模型
G1d =0.004555 z + 0.00178----------------------z^2 - 1.05 z + 0.04979% G2的离散模型
G2d =0.04117 z + 0.03372----------------------z^2 - 1.511 z + 0.5488% 闭环离散模型
GHd =0.0001875 z^2 + 0.0002268 z + 6e-05------------------------------------------------z^4 - 2.561 z^3 + 2.185 z^2 - 0.6512 z + 0.02738

3.2.3 实战3

实验要求：使用$\mathrm{MATLAB}$绘制离散系统$G(z)=\frac{2z^2-3.4z+1.5}{z^2-1.6z+0.8}$的带栅格线的根轨迹图。

解：

% 实例Chapter11.3.2.3
clc;clear;num=[2,-3.4,1.5];den=[1,-1.6,0.8];
zgrid('new');rlocus(num,den);
set(findobj(get(gca,'Children'),'LineWidth',0.5),'LineWidth',1.5);
title('离散系统的根轨迹图','FontSize',15);

3.2.4 实战4

实验要求：已知系统闭环传递函数为：$G(z)=\frac{z^3-1.3z^2+1.22z+0.51}{z^4+0.522z^3+0.4z^2+0.0086z-0.3915}，T=0.5$，绘制此系统的零极点图，判断此系统的稳定性。

解：

% 实例Chapter11.3.2.4
clc;clear;% 建立系统模型
num=[1,-1.3,1.22,0.51];
den=[1,0.522,0.4,0.0086,-0.3915];
G=tf(num,den);% 绘制零极点图
pzmap(G);
axis([-1.2,1.2,-1.2,1.2]);axis equal;
set(findobj(get(gca,'Children'),'LineWidth',0.5),'LineWidth',1.5);
title('离散系统的零极点图','FontSize',15);

• 由零极点图可知，闭环传递函数所有极点都位于单位圆内部，闭环系统稳定；

3.2.5 实战5

实验要求：若某控制系统结构如下图所示，其中：$D_1(z)=\frac{3.4z^{-1}-1.5z^{-2}}{1-1.6z^{-1}+0.8z^{-2}}$，$G_1$是零阶保持器，$G_1(s)=\frac{1-{\mathrm{e}}^{-0.05s}}{s}，G_2(s)=\frac{0.25}{s^2+3s+2}$，采样周期：$T_s=0.05\mathrm{s}$，求系统开环和闭环的$z$传递函数，及$s$传递函数，当输入为单位阶跃函数时，求其输出。

解：

【求系统开环和闭环$s$传递函数】

% 实例Chapter11.3.2.5
clc;clear;% 系统开环和闭环s传递函数
dnum1=[3.4,-1.5];dden1=[1,-1.6,0.8];Ts=0.05;
sysd1=tf(dnum1,dden1,Ts);           % D1的z传递函数模型
sysc1=d2c(sysd1,'zoh');num2=[0.25];den2=[1,3,2];sys2=tf(num2,den2);
sysc2=sysc1*sys2                    % 开环传递函数
sysbc=feedback(sysc2,1)             % 闭环传递函数[num,den]=tfdata(sysbc,'v');        % 提取闭环传递函数的分子和分母
p=roots(den)

% 结果显示：% 开环s传递函数
sysc2 =13.99 s + 216---------------------------------------------s^4 + 7.463 s^3 + 106.4 s^2 + 281.8 s + 181.9% 闭环s传递函数
sysbc =13.99 s + 216-------------------------------------------s^4 + 7.463 s^3 + 106.4 s^2 + 295.8 s + 398% 特征方程的根
% 由于特征方程的根均在s左半平面，因此，系统是稳定的。
p =-2.1667 + 9.1429i-2.1667 - 9.1429i-1.5647 + 1.4351i-1.5647 - 1.4351i

【求取系统开环和闭环$z$传递函数】

% 系统开环和闭环z传递函数
dnum1=[3.4,-1.5];dden1=[1,-1.6,0.8];Ts=0.05;
sysd1=tf(dnum1,dden1,Ts);           % D1的z传递函数模型num2=[0.25];den2=[1,3,2];
sys2=tf(num2,den2);
sysd2=c2d(sys2,Ts,'zoh');           % G1和G2串联的z传递函数模型sysd=sysd1*sysd2                    % 系统开环z传递函数
sysbd=feedback(sysd,1)              % 系统闭环z传递函数[dnum,dden]=tfdata(sysbd,'v');
pd=roots(dden)                      % 闭环系统特征根t=0:0.05:5;
y=dstep(dnum,dden,101);
stem(t,y);                          % 棒图显示响应曲线
xlabel('t');ylabel('y');
grid;
set(findobj(get(gca,'Children'),'LineWidth',0.5),'LineWidth',1.5);
title('控制系统响应曲线','FontSize',15);

% 结果显示：% 开环z传递函数
sysd =0.001011 z^2 + 0.0005156 z - 0.0004242---------------------------------------------z^4 - 3.456 z^3 + 4.63 z^2 - 2.862 z + 0.6886% 闭环z传递函数
sysbd =0.001011 z^2 + 0.0005156 z - 0.0004242----------------------------------------------z^4 - 3.456 z^3 + 4.631 z^2 - 2.861 z + 0.6881% 特征方程的根
pd =0.8030 + 0.3935i0.8030 - 0.3935i0.9250 + 0.0697i0.9250 - 0.0697i

3.2.6 实战6

实验要求：某控制系统结构如下图所示，其中：$D_1(z)=\frac{3.4z^{-1}-1.5z^{-2}}{1-1.6z^{-1}+0.8z^{-2}}$，$G_1$是零阶保持器，$G_1(s)=\frac{1-{\mathrm{e}}^{-0.05s}}{s}，G_2(s)=\frac{0.25}{s^2+3s+2}$，采样周期：$T_s=0.05\mathrm{s}$，求系统频率特性参数，绘制系统的$\mathrm{Bode}$图、$\mathrm{Nyquist}$图和$\mathrm{Nichols}$图。

解：

% 实例Chapter11.3.2.6
clc;clear;% 系统开环和闭环s传递函数
dnum1=[3.4,-1.5];dden1=[1,-1.6,0.8];Ts=0.05;
sysd1=tf(dnum1,dden1,Ts);           % D1的z传递函数模型num2=[0.25];den2=[1,3,2];sys2=tf(num2,den2);
sysd2=c2d(sys2,Ts,'zoh');
sysd=sysd1*sysd2                    % 开环z传递函数[dnumc,ddenc]=tfdata(sysd,'v');     % 提取开环传递函数的零极点figure(1);
[Gm,Pm,Wcg,Wcp]=margin(sysd)        % 求系统频率特性参数
margin(sysd);grid;
set(findobj(get(gca,'Children'),'LineWidth',0.5),'LineWidth',1.5);
title('控制系统margin图','FontSize',15);w=0.01:0.01:100;figure(2);
dnyquist(dnumc,ddenc,Ts,w);grid;
set(findobj(get(gca,'Children'),'LineWidth',0.5),'LineWidth',1.5);
title('控制系统nyquist图','FontSize',15);figure(3);
nichols(sysd);grid;
set(findobj(get(gca,'Children'),'LineWidth',0.5),'LineWidth',1.5);
title('控制系统nichols图','FontSize',15);

% 结果显示：% 控制系统开环模型
sysd =0.001011 z^2 + 0.0005156 z - 0.0004242---------------------------------------------z^4 - 3.456 z^3 + 4.63 z^2 - 2.862 z + 0.6886% 频率特性参数
Gm =10.8194Pm =133.7794Wcg =6.4303Wcp =0.5625

【控制系统$\mathrm{Bode}$图】

【控制系统$\mathrm{Nyquist}$图】

【控制系统$\mathrm{Nichols}$图】

3.3 综合实例

实验要求：给定单位负反馈离散控制系统，其采样周期为$1\mathrm{s}$，开环传递函数为：$G(s)=\frac{s+1}{s^2}$与零阶保持器$\mathrm{ZOH}$串联，开环增益为$K$。求闭环控制系统稳定的条件，且绘制$K$取不同值时闭环系统的阶跃响应曲线。

解：

【$STEP1$：建立控制系统的数学模型】

% 实例Chapter11.3.3
clc;clear;% 建立控制系统模型
Ts=1;sys_K=1;
num=[1,1];den=[1,0,0];
sysc=tf(num,den);                   % 连续系统传递函数
sysd=c2d(sysc,Ts,'zoh');            % 离散系统传递函数
sys_open=sys_K*sysd                 % 系统的开环传递函数

% 系统开环传递函数
sys_open =1.5 z - 0.5-------------z^2 - 2 z + 1

【$STEP2$：绘制系统的根轨迹】

% 绘制系统的根轨迹
figure(1);
rlocus(sysd);axis([-2,2,-1,1]);
set(findobj(get(gca,'Children'),'LineWidth',0.5),'LineWidth',1.5);
title('控制系统根轨迹图','FontSize',15);% 验证临界稳定的K值
% K=2时是临界稳定值
sys_K=2;
figure(2);
margin(sys_K*sysd);
set(findobj(get(gca,'Children'),'LineWidth',0.5),'LineWidth',1.5);
title('控制系统margin图','FontSize',15);figure(3);
[dnum,dden]=tfdata(sys_K*sysd,'v');
dnyquist(dnum,dden,Ts);axis([-5,5,-2.5,2.5]);
grid on;
set(findobj(get(gca,'Children'),'LineWidth',0.5),'LineWidth',1.5);
title('控制系统Nyquist图','FontSize',15);

【控制系统根轨迹图】

【控制系统$\mathrm{Bode}$图】

【控制系统$\mathrm{Nyquist}$图】

【$STEP3$：控制系统阶跃响应】

% 绘制不同K值的阶跃响应
sys_K=1;
figure(4);
sys_close=feedback(sys_K*sysd,1);
[dnumc,ddenc]=tfdata(sys_close,'v');
dstep(dnumc,ddenc,25);
set(findobj(get(gca,'Children'),'LineWidth',0.5),'LineWidth',1.5);
title('K=1时控制系统阶跃响应','FontSize',15);sys_K=2;
figure(5);
sys_close=feedback(sys_K*sysd,1);
[dnumc,ddenc]=tfdata(sys_close,'v');
dstep(dnumc,ddenc,25);
set(findobj(get(gca,'Children'),'LineWidth',0.5),'LineWidth',1.5);
title('K=2时控制系统阶跃响应','FontSize',15);sys_K=3;
figure(6);
sys_close=feedback(sys_K*sysd,1);
[dnumc,ddenc]=tfdata(sys_close,'v');
dstep(dnumc,ddenc,25);
set(findobj(get(gca,'Children'),'LineWidth',0.5),'LineWidth',1.5);
title('K=3时控制系统阶跃响应','FontSize',15);

• 当$K=1$时，闭环系统稳定，阶跃响应曲线收敛；
• 当$K=2$时，闭环系统临界稳定，阶跃响应曲线等幅振荡；
• 当$K=3$时，闭环系统不稳定，阶跃响应曲线发散；

【完整代码】

% 实例Chapter11.3.3
clc;clear;% 建立控制系统模型
Ts=1;sys_K=1;
num=[1,1];den=[1,0,0];
sysc=tf(num,den);                   % 连续系统传递函数
sysd=c2d(sysc,Ts,'zoh');            % 离散系统传递函数
sys_open=sys_K*sysd                 % 系统的开环传递函数% 绘制系统的根轨迹
figure(1);
rlocus(sysd);axis([-2,2,-1,1]);
set(findobj(get(gca,'Children'),'LineWidth',0.5),'LineWidth',1.5);
title('控制系统根轨迹图','FontSize',15);% 验证临界稳定的K值
sys_K=2;
figure(2);
margin(sys_K*sysd);
set(findobj(get(gca,'Children'),'LineWidth',0.5),'LineWidth',1.5);
title('控制系统margin图','FontSize',15);figure(3);
[dnum,dden]=tfdata(sys_K*sysd,'v');
dnyquist(dnum,dden,Ts);axis([-5,5,-2.5,2.5]);
grid on;
set(findobj(get(gca,'Children'),'LineWidth',0.5),'LineWidth',1.5);
title('控制系统Nyquist图','FontSize',15);% 绘制不同K值的阶跃响应
sys_K=1;
figure(4);
sys_close=feedback(sys_K*sysd,1);
[dnumc,ddenc]=tfdata(sys_close,'v');
dstep(dnumc,ddenc,25);
set(findobj(get(gca,'Children'),'LineWidth',0.5),'LineWidth',1.5);
title('K=1时控制系统阶跃响应','FontSize',15);sys_K=2;
figure(5);
sys_close=feedback(sys_K*sysd,1);
[dnumc,ddenc]=tfdata(sys_close,'v');
dstep(dnumc,ddenc,25);
set(findobj(get(gca,'Children'),'LineWidth',0.5),'LineWidth',1.5);
title('K=2时控制系统阶跃响应','FontSize',15);sys_K=3;
figure(6);
sys_close=feedback(sys_K*sysd,1);
[dnumc,ddenc]=tfdata(sys_close,'v');
dstep(dnumc,ddenc,25);
set(findobj(get(gca,'Children'),'LineWidth',0.5),'LineWidth',1.5);
title('K=3时控制系统阶跃响应','FontSize',15);

$$=========================\mathrm{COMPLETE}========================$$