《MATLAB/Simulink系统仿真超级学习手册》——2.7 MATLAB编程实例

本节书摘来自异步社区《MATLAB/Simulink系统仿真超级学习手册》一书中的第2章，第2.7节，作者：MATLAB技术联盟 , 石良臣著，更多章节内容可以访问云栖社区“异步社区”公众号查看

2.7 MATLAB编程实例

MATLAB/Simulink系统仿真超级学习手册
MATLAB强大的运算功能和图形功能，是实现各领域的编程仿真的有力工具。本节将综合本章所讲的内容介绍MATLAB的一些编程实例，使读者对MATLAB语言的应用有一些基本的了解。

2.7.1 汉诺塔问题
Hanoi（汉诺）塔问题，是一个古典的数学问题，是一个用递归方法解题的典型例子。问题是这样的：古代有一个梵塔，塔内有3个座A、B、C，开始时A座上有n个盘子，盘子大小不等，大的在上，小的在下如图2-35所示。

现在，要将这n个盘子从A座移到C座，但每次只允许移动一个盘子，且在移动过程中在3个座上都保持大盘在下，小盘在上。在移动过程中可以利用B座。

显然，读者是不大可能直接写出移动盘子的每一个具体步骤的，这里给出该问题的递归解法：

（1）将A上n-1个盘子借助C座先移动到B座上；

（2）把A座上剩下的一个盘子移动到C座上；

（3）将n-1个盘子从B座借助于A座移动到C座上。

上面第（1）步和第（3）步，都是把n-1个盘子从一个座移到另一个座上，采取的办法是一样的，只是座的名字不同而已。

为使之一般化，可以将第（1）步和第（3）步表示为：将“one”座上n-1个盘子移动到“two”座（借助“three”座）。只是在第（1）步和第（3）步中，one、two、three和A、B、C的对应关系不同。

对于第（1）步，对应关系是，one对应A，two对应B，three对应C。对于第（3）步，对应关系是，one对应B，two对应C，three对应A。

因此，可以把上面3个步骤分成两类操作：

（1）将n-1个盘子从一个座移动到另一个座上（n>1）；

（2）将1个盘子从一个座上移动到另一个座上。

分别用两个函数实现上面的两类操作，用hanoi函数实现上面第一类操作，用move函数实现上面第二类操作。

函数调用hanoi(n,one,two,three)表示将n个盘子从“one”座移动到“three”座的过程（借助“two”座），函数调用move(x,y)表示将1个盘子从x座移动到y座的过程。x和y是代表A、B、C座之一，根据每次不同情况分别取A、B、C代入。

编写函数文件，输入下面的语句，保存为epHanoi.m。

function epHanoi
% 汉诺塔递归算法
% n是需要移动的盘子个数
% 最后打印出全部移动步骤
n=input('请输入盘子个数：');
disp(['移动',num2str(n),'个盘子的步骤为：'])
(n,'A','B','C')
function (n,one,two,three)
if n==1
　　move(one,three);
else
　　hanoi(n-1,one,three,two);
　　move(one,three);
　　hanoi(n-1,two,one,three);
end
function move(x,y)
disp([x,'->',y])

在MATLAB命令窗口输入epHanoi命令，调用函数epHanoi，程序会先提示“请输入盘子个数：”，读者可以任意输入一个整数（如3），按回车键确认后，出现下面的结果。

>> epHanoi
请输入盘子个数：3
移动3个盘子的步骤为：
A->C
A->B
C->B
A->C
B->A
B->C
A->C

在这个例子中move函数并未真正移动盘子，而只是输出了移动盘子的方案。

注意：

（1）函数文件可以没有输入，也可以没有输出；

（2）从这个例子可以看出，MATLAB支持递归运算；

（3）函数文件epHanoi.m的代码中，出现了三个function，第一个为主函数，后两个为子函数，外部函数无法访问这两个子函数，除非将其定义为global全局的。一个函数下的子函数可以互相调用。

MATLAB在自动控制中的应用
【例2-28】已知单位负反馈系统的开环传递函数为

画出K=0~∞时的闭环系统根轨迹，并求出临界时的K值及闭环极点。

函数rlocus可以计算单输入单输出开环系统根的位置，进而给出其单位负反馈系统的闭环根轨迹，其基本调用格式为：

rlocus(sys)

其中，sys表示开环系统，一般用函数tf(num,den)来创建。num、den分别表示传递函数的分子、分母，用对应多项式的系数构成的向量表示。

编写脚本文件，输入下面的语句，保存为ep2_28.m。

d1=[0.05 1];
d2=[0.05 0.2 1];
den1=conv(d1,d2);
den=[den1 0];
num=1;
sys=tf(num,den);
rlocus(sys)
rlocfind(sys)

运行脚本文件ep2_28.m。程序执行到rlocus(sys)时，将出现如图2-36所示的根轨迹图。

注意：conv(d1,d2)是求向量d1和向量d2的卷积。

为了计算K临界值和此时的闭环极点，在根轨迹绘制完成（执行函数rlocus）后，再执行rlocfind命令。

程序执行rlocfind(sys)语句后，图形窗口Figure 1上将产生一个跟随鼠标移动的十字形指针，选中轨迹与虚轴的焦点位置单击鼠标左键，MATLAB命令窗口中就会显示出闭环根以及参数K的值。例2-28中，所求临界值K=3.49。

函数roots可以用于求多项式的根，在MATLAB命令窗口输入roots([den1,3.49])，求得闭环极点为s1=-19.7878、s2=-4.2244、s3=0.0061+4.0866j和s4=0.0061-4.0866j。

MATLAB在电力信号分析处理中的应用
【例2-29】利用函数生成一组数据用以模拟电力现场测量数据，并利用MATLAB编程实现其频谱分析。其中

在信号x中加入随机信号，用以模拟信号噪声。

编写脚本文件，输入下面的语句，保存为ep2_29.m。

clear
fs=500;t=0:1/fs:0.2;
f1=50;f2=200;
x=2*sin(2*pi*f1*t)+cos(2*pi*f2*t);
subplot(4,1,1);plot(x)
title('f1（50Hz）\f2（200Hz）的正弦信号，初相0')
xlabel('序列（n）')
grid on
number=512;
y=fft(x,number);
n=0:length(y)-1;
f=fs*n/length(y);
subplot(4,1,2);plot(f,abs(y))
title('f1\f2的正弦信号的FFT（512点）')
xlabel('频率Hz')
grid on
x=x+randn(1,length(x));
subplot(4,1,3);plot(x)
title('f1\f2的正弦信号（含随机噪声）')
xlabel('序列（n）')
grid on
y=fft(x,number);
n=0:length(y)-1;
f=fs*n/length(y);
subplot(4,1,4);plot(f,abs(y))
title('f1\f2的正弦信号（含随机噪声）的FFT（512点）')
xlabel('频率Hz')
grid on

运行脚本文件ep2_29.m，运行结果如图2-37所示。

注意：

（1）在脚本文件开头处加入clear命令，是为了清除之前存留在工作空间中的变量和数据，避免与之后新定义或要使用的变量发生混乱；

（2）函数randn可以产生一个服从标准分布的随机矩阵。