第十二章 MATLAB的学科应用

第12章MATLAB的学科应用【本章学习目标】●综合运用所学MATLAB程序设计知识，学会针对具体的问题，选择合适的解题方案，灵活掌握利用MATLAB解决实际问题的方法。●了解MATLAB中几个典型的学科工具箱。●熟悉MATLAB在相关学科领域的应用案例。MATLAB有40多个工具箱，大致可分为两类：功能性工具箱:主要用来扩充MATLAB的符号计算功能、图形建模仿真功能、文字处理功能以及与硬件的实时交互功能，能用于多种学科。学科性工具箱:涵盖了控制系统设计与分析、数字信号处理、数字图像处理、金融财务分析、生物遗传过程等许多专业领域，可以利用这些工具箱进行相关领域的科学研究，解决相关领域的实际问题。12.1MATLAB在优化问题中的应用在日常生活和实际工作中，人们对于同一个问题往往会提出多个解决方案，并通过各方面的论证从中提取最佳方案。最优化方法就是专门研究如何从多个方案中科学合理地提取出最佳方案的科学。最优化方法的应用和研究已经深入到了生产和科研的各个领域，如土木工程、机械工程、化学工程、运输调度、生产控制、经济规划、经济管理等，并取得了显著的经济效益和社会效益。12.1.1优化模型与优化工具1．优化模型优化模型是用数学关系式表示的研究对象的某种本质特征，它有如下3个要素:●决策变量：问题中要确定的未知量，用以表示优化方案。●目标函数：表示所待决策问题期望达到的目标，它是决策变量的函数。●约束条件：指决策变量取值时受到的各种资源条件的限制，通常用含决策变量的等式或不等式表示。2．优化函数

MATLAB的优化工具箱（OptimizationToolbox）提供了一组优化函数，在使用这组优化函数时，用户定义的待优化函数作为输入参数传递给这些函数。3．优化工具图形用户界面使用优化工具图形用户界面，可以通过鼠标操作选择优化参数和运行优化。在MATLAB的命令窗口输入“optimtool”命令可以打开优化工具图形用户界面12.1.2应用实例利用MATLAB的优化工具箱，可以求解线性规划、非线性规划和多目标规划问题。使用优化工具箱，先要定义目标函数和约束条件，然后设置优化参数（如算法），最后调用优化工具求解。1．一元函数最小值问题fminbnd函数用于求解有约束的一元函数最小值问题，其调用格式为[x,fval]=fminbnd(fun,x1,x2,options)其中，fval为目标函数的最小值，fun为目标函数，解的约束为x1≤x≤x2。【例12.1】对边长为3m的正方形铁板，在4个角剪去相等的正方形以制成方形无盖水槽，问如何剪法使水槽的容积最大？设剪去的正方形的边长为x，则水槽的容积为(3

−

2x)2x。fminbnd函数用于求最小值，因此将求解水槽最大容积转换为miny

=

−

(3

−

2x)2x，0

<

x

<

1.5在MATLAB命令窗口输入命令：[x,fval]=fminbnd(@(x)-(3-2*x)^2*x,0,1.5);求解得到：x=0.5000fval=−2.0000即剪掉的正方形的边长为0.5m时，水槽的容积最大，最大容积为2m3。此题也可使用优化工具的GUI来求解。【例12.2】生产计划问题。某企业在计划期内计划生产甲、乙、丙3种产品。这些产品分别需要在设备A、B上加工，需要消耗材料C、D，按工艺资料规定，单件产品在不同设备上加工及所需要的资源如表12.2所示。已知在计划期内设备的加工能力各为200台时，可供材料分别为360kg、300kg；每生产一件甲、乙、丙3种产品，企业可获得利润分别为40、30、50元，假定市场需求无限制。企业决策者应如何安排生产计划，使企业在计划期内总的利润收入最大？设在计划期内生产这3种产品的产量为x1、x2、x3，用Z表示利润，则有Z

=

40x1

+

30x2

+

50x3。在安排3种产品的计划时，不得超过设备A、B的可用工时，材料消耗总量不得超过材料C、D的供应量，生产的产量不能小于零。企业的目标是要使利润达到最大，这个问题的数学模型为表12.2单位产品资源消耗产品资源甲乙丙现有资源设备A312200设备B224200材料C451360材料D235300利润（元/件）4030503．非线性规划求解非线性规划问题的函数是fmincon，其调用格式为[x,fval]=fmincon(f,x0,A,b,Aeq,beq,lb,ub)其中，x0是初值，其余参数含义与linprog函数的相同。【例12.3】设有400万元资金，要求4年内使用完，若在一年内使用资金x万元，则可得效益万元（效益不能再使用），当年不用的资金可存入银行，年利率为10%。试制定出资金的使用计划，以使4年效益之和为最大。

设变量xi表示第i年所使用的资金数，则有（1）定义目标函数functionf=xymb(x)f=-(sqrt(x(1))+sqrt(x(2))+sqrt(x(3))+sqrt(x(4)));（2）主程序youh3.m为x0=[1;1;1;1];A=[1000;1.1100;1.211.110;1.3311.211.11];b=[400;440;484;532.4];[x,fval]=fmincon(@xymb,x0,A,b)*12.2MATLAB在控制系统中的应用控制系统的分析工具包括控制系统工具箱（ControlSystemToolbox）、系统辨识工具箱（SystemIdentificationToolbox）、模糊逻辑工具箱（FuzzyLogicToolbox）、鲁棒控制工具箱（RobustControlToolbox）、模型预测控制工具箱（ModelPredictiveControlToolbox）、线性矩阵不等式工具箱（LMIControlToolbox）等。12.2.1控制系统工具箱控制系统工具箱包含了丰富的线性系统分析和设计函数，并以LTI对象为基本数据类型对线性时不变系统进行操作与控制。控制系统工具箱能够完成系统的时域和频域分析。1．系统分析控制系统的分析包括系统的时域分析、频域分析、稳定性分析及根轨迹分析。命令行方式调用各种函数完成模型分析。2．系统设计控制系统工具箱支持4种线性模型表述方式：传递函数模型、零极点增益模型、状态空间模型和频域响应数据模型。每一种模型形式都有相应的LTI对象的表达。除了模型数据以外，LTI对象还可以存储离散系统的采样时间、时间延迟、输入输出名称、模型注释等。控制系统工具箱提供了GUI和命令行两种方式进行控制系统的设计。命令行方式调用各种函数建立系统模型SISO（Single-InputSingle-Output）设计器是控制系统工具箱所提供的线性系统设计器，它为用户设计单输入单输出线性控制系统提供了非常友好的图形界面。在SISO设计器中，用户可以同时使用根轨迹图与波特图，通过修改线性系统零点、极点以及增益等传统设计方法进行SISO线性系统设计。在MATLAB命令窗口中输入sisotool命令启动SISO设计器。使用控制系统工具箱中的函数还可以进行各种系统的补偿设计，如LQG（线性二次型设计）、RootLocus（线性系统的根轨迹设计）、Poleplacement（线性系统的极点配置）、Observer-basedregulator（线性系统观测器设计）等12.2.2应用实例G1=tf([251],[123]);G2=zpk([0-2],-10,1);Cloop=feedback(G1,G2); %构造闭环系统的传递函数step(Cloop); %求系统的阶跃响应holdon;impulse(Cloop,'r-.'); %求系统的冲激响应legend('Step','Impulse','--','r-.');【例12.5】在图12.6所示的控制系统中，Gc(s)是一比例环节。试利用MATLAB控制系统工具箱求取系统临界稳定的开环增益值K。

（1）建立系统模型，然后启动SISO设计器。在MATLAB命令窗口输入以下命令：Gp=tf(1,[17100]);sisotoolSISO设计器通过零极点配置、根轨迹分析、系统波特图分析等传统的方法对线性系统进行设计。（2）导入系统数据。（3）求临界稳定增益。*12.3MATLAB在信号处理中的应用12.3.1信号处理工具箱与Simulink模块集在MATLAB中，与信号处理有关的组件包括信号处理工具箱（SignalProcessingToolbox）、神经网络工具箱（NeuralNetworkToolbox）、小波分析工具箱（WaveletToolbox）、统计工具箱（StatisticsToolbox）、通信工具箱（CommunicationsToolbox）、图像处理工具箱（ImageProcessingToolbox）等以及Simulink中的各种信号处理模块。1．信号处理工具箱信号处理工具箱是一个信号处理的工业级标准算法集合，利用它可以完成数字或模拟系统中常规的信号处理任务。例如，信号生成、信号重采样、各种积分变换、模拟/数字滤波器的设计和分析、频谱分析、随机信号分析、参数模型分析、线性预测等。信号处理工具箱提供的函数包括以下几类。（1）产生基本信号的函数。例如，方波函数square、锯齿波函数sawtooth、矩形脉冲函数tectpuls、三角脉冲函数tripuls、高斯调制正弦波脉冲函数gauspuls、扫频余弦信号函数chrip等。这些基本信号是信号处理的基础。（2）信号变换函数。例如，快速傅立叶变换函数fft、离散余弦变换函数dct、Hilbert变换函数hilbert和ChripZ变换函数czt等。（3）数字滤波器的设计函数。提供了一套完整的有限冲激响应（FIR）和无限冲激响应（IIR）方法，这些方法支持低通、高通、带通、带阻和多频带滤波器的设计与分析。（4）随机信号的分析函数。随机信号的分析包括信号的相关性分析和功率谱估计。FDATool（FilterDesignandAnalysis）是设计常规滤波器的图形用户界面，它涵盖了信号处理工具箱中所有的滤波器设计方法，利用它可以分别设计出满足新年指标的滤波器，并可查看该滤波器的各种分析图形。在MATLAB命令窗口输入fdatool命令打开FDATool界面。SPTool（SignalProcessingTool）是一个信号处理的GUI集成环境，提供对信号、滤波器和频谱分析函数的访问入口。在MATLAB命令窗口输入“sptool”命令，打开SPTool图形用户界面。利用SPTool可以从MATLAB工作空间或从文件中导入已经设计好的信号和滤波器频谱进行时域和频域分析，还可以用于设计滤波器和实现对信号的滤波，以及对随机信号进行谱分析。2．Simulink模块集Simulink的信号处理模块集（SignalProcessingBlockset）也提供了一系列信号处理模块，允许用户在不进行底层编程的情况下设计和仿真实时系统。信号处理模块约定以列向量表示单通道信号，在多通道情况下，每一列代表一个通道，每一行对应一个采样点。12.3.2应用实例【例12.6】采用海明窗函数设计一个95阶的线性相位带通滤波器，其技术指标要求为：采样频率1

000Hz，通带下限截止频率为70Hz，通带上限截止频率为84Hz设计。n=95;fc1=70;fc2=84;fs=1000;w1=2*pi*fc1/fs;w2=2*pi*fc2/fs;window=hamming(n+1);h=fir1(n,[w1/piw2/pi],window);freqz(h,1,512);【例12.7】已知信号x(t)

=

sin(2

f1t)

+

2sin(2

f2t)

+

w(t)，f1

=

50Hz，f2

=

120Hz，W(t)为标准差是0.1的白噪声，采样频率为1

000Hz，信号长度为1

024。使用周期图法求信号的功率谱。周期图法是将信号的采样数据进行Fourier变换求取功率谱密度估计的方法。利用MATLAB信号处理工具箱中已经定义的周期图法函数，编写如下程序：f=[50;120];A=[12];fs=1000;N=1024;n=0:N-1;t=n/fs;x=A*sin(2*pi*f*t)+0.1*randn(1,N);%生成带噪声的信号periodogram(x,[],N,fs);%绘制功率谱Pxx=periodogram(x,[],'twosided',N,fs);Pow=(fs/length(Pxx))*sum(Pxx);%计算平均功率12.4MATLAB在数字图像处理中的应用图像处理工具包括图像处理工具箱（ImageProcessingToolbox）、图像获取工具箱（ImageAcquisitionToolbox）、地理信息处理工具箱（MappingToolbox）、信号处理工具箱、统计工具箱等，提供图像数据采集、图像处理和分析的工具。12.4.1图像处理工具箱MATLAB图像处理工具箱支持4种图像类型，分别为真彩色图像、索引色图像、灰度图像和二值图像。1．常用图像操作函数常用的图像操作函数包括读取图像的函数imread、输出图像的函数imwrite、在屏幕上显示图像的函数imshowimage等，对图像进行裁剪可用imcrop函数，图像的插值缩放可用imresize函数实现，旋转用imrotate函数实现。2．图像增强功能图像增强采用一系列技术去改善图像的视觉效果或将图像转换成一种更适合于人眼观察和机器自动分析的形式。常用的图像增强方法有以下几种。（1）灰度直方图均衡化。histeq函数（2）灰度变换法。imadjust函数（3）平滑与锐化滤波。用fspecial函数创建预定义的滤波算子，然后用filter2或conv2函数在实现卷积运算的基础上进行滤波。3．边缘检测和图像分割功能MATLAB工具箱提供的edge函数可针对sobel算子、prewitt

算子、Roberts算子、log算子和canny算子实现检测边缘的功能。4．图像变换功能图像处理工具箱提供了图像变换函数，如fft2与ifft2函数分别实现二维快速傅立叶变换与逆变换，dct2与idct2函数实现二维离散余弦变换与逆变换，wavedec2和waverec2函数分别实现二维多尺度小波变换与逆变换。5．图像获取12.4.2应用实例【例12.8】显示两幅图像在空域和频域融合后的效果。图12.16原图像图12.17融合图像X1=imread(‘pic_1.jpg');X2=imread(‘pic_2.jpg');subplot(1,2,1);image(X1);title('图像1')subplot(1,2,2);image(X2);title('图像2')figuresubplot(1,3,1);image((X1+X2)/2);%在空域内直接融合title('直接相加融合')F1=fft2(X1);F2=fft2(X2);XR1=uint8(abs(ifft2(F1+F2)/2));subplot(1,3,2);image(XR1);%显示经过FFT变换融合后的图像title('傅立叶变换融合')%分别对两幅原图像进行小波分解[C1,L1]=wavedec2(X1,2,'sym4');[C2,L2]=wavedec2(X2,2,'sym4');C=C1+C2;%对分解系数进行融合XR2=waverec2(C,L1,'sym4');%对融合后的信号进行图像重构subplot(1,3,3)image(uint8(XR2/2));%显示经过小波变换融合后的图像title('小波变换融合')12.5MATLAB在经济和金融领域中的应用经济和金融方面的建模和分析工具提供财务数据分析、投资组合最佳化、商品评价、风险管理、敏感度分析，以及经济模型建立与预测的工具。12.5.1经济和金融领域的工具箱经济工具箱金融工具箱固定收益工具箱统计工具箱财务数据来源工具箱神经网络工具箱优化工具箱12.5.2应用实例【例12.9】假设投资者有两种总资产，总价值为1亿元，资产权重分别为30%与70%，资产的日波动率均值分别为0.002和0.004，标准差为0.03和0.01。这两种资产的相关系数为0.8，给定置信度为0.99，求该资产在30天的投资组合风险值。w=[0.3,0.4];ret=[0.002,0.004];d=[0.03,0.01];cov=[1,0.8;0.8,1];time=30;pret=time*dot(w,ret);%计算总收益sd=w.*d;pd=sqrt(sd*cov*sd'*time);%计算总资产的标准差Var=portvrisk(pret,pd,0.01,1);程序执行后，得到Var=0.0

924。由此即可求出该资产的投资组合风险值为0.0

924亿元。【例12.10】某商店为了确定向公司A或公司B购买某种产品，将A、B公司以往各次进货的次品率进行比较，数据如表12.4所示，设两样本独立。问两公司商品的质量有无显著差异。设两公司的商品的次品密度最多只差一个平移，取α=

0.05。a=[7.0,3.5,9.6,8.1,6.2,5.1,10.4,4.0,2.0,10.5];b=[5.7,3.2,4.2,11.0,9.7,6.9,3.6,4.8,5.6,8.4,10.1,5.5,12.3];[p,h]=ranksum(a,b)p=0.8282h=0运行结果表明两样本总体均值相等的概率为0.8282，并不很接近于零，且h

=

0说明可以接受原假设，即认为两个公司的商品的质量无明显差异。补充实例1：实验数据处理与作图x1=[200350500650800950110012501400155017001900];x2=[265026501842184210341034];x3=[280026003000270025002900260030002800250029002700];x4=[875497649865];y=[0.1510.1130.1990.1160.0910.1420.0990.1350.1280.0290.1160.016]';>>X=[ones(12,1)x1'x2'x3'x4'(x1.^2)'(x2.^2)'(x3.^2)'(x4.^2)'];[b,bint,r,rint,stats]=regress(y,X);>>b,statsb=0.383603491362780.00001000504712-0.00332404632897-0.000352943401230.01421001562287-0.000000035841020.000040342402350.00000009851594-0.00107630787579stats=0.9944025921256666.620296469186770.00272063036105将计算得的b阵代入可算得最优条件为TC=140℃、tc=41s、Ta=1790℃，ta=7s补充实例2：药学建模例3药代动力学实验数据处理及参数估计%输入数据

t=[251020306090120150180210];

c=[766659470300239180160145132123112];%x0为给A,alpha,B,beta初值

x0=[0000];%定义模型函数f=inline('x(1)*exp(-x(2)*t)+x(3)*exp(-x(4)*t)','x','t');%进行非线性曲线拟合

[b,r,j]=nlinfit(t,c,f,x0);b=b';%得到各个系数值

A=b(1),alpha=b(2),B=b(3),beta=b(4)记中心室和吸收室的容积分别为V,V1，而t时刻的血药浓度分别为c(t),c1(t);中心室的排除速率为k，吸收速率为k1(这里k和k1分别是中心室和吸收室血药浓度变化率与浓度本身的比例系数，设t=0时刻口服剂量为d的药物，容易写出吸收室的血药浓度c1(t)的微分方程为

(dc_1)/dt=-k_1*c_1

c_1(0)=d/V_1中心室血药浓度c(t)的变化率由两部分组成：与c成正比的排除（比例系数k）;与c1成正比的吸收（比例系数k1）.再考虑到中心室和吸收室的容积分别为V,V1,得到c(t)的微分方程为

dc/dt=-k*c+V_1/V*k_1*c_1

c(0)=0由以上两个微分方程不难结出中心室血药浓度

c(t)=(d/V)*(k_1/(k_1-k))*(e^(-k*t)-e^(-k_1*t)).在制定给药方案时必须知道这种药物的3个参数k1,k,b(=d/V),实际中通常通过实验数据确定。设t=0时刻口服一定剂量的药物，表7.11是实验数据c(t),请由此确定k,k1,b.t0.0830.1670.250.500.751.01.52.253.04.06.08.010.012.0c(t)10.921.127.336.435.538.434.824.223.615.78.28.32.21.8模型与解答：这是一个非线性拟合问题，可以用lsqcurvefit求解。MATLAB程序如下：

functionf=ex08(x,t)

f=x(3)*x(1)/(x(1)-x(2))*(exp(-x(2)*t)-exp(-x(1)*t));x0=[2,3,4];

t=[0.083,0.167,0.25,0.50,0.75,1.0,1.5,2.25,3.0,4.0,6.0,8.0,10.0,12.0];

c=[10.9,21.1,27.3,36.4,35.5,38.4,34.8,24.2,23.6,15.7,8.2,8.3,2.2,1.8];

[x1,norm1,res]=lsqcurvefit(@ex08,x0,t,c);

[x2,norm2,res]=lsqcurvefit(@ex08,x1,t,c);

[x3,norm3,res]=lsqcurvefit(@ex08,x2,t,c);

[x4,norm4,res]=lsqcurvefit(@ex08,x3,t,c);

[x5,norm5,res]=lsqcurvefit(@ex08,x4,t,c);[x20,norm20,res]=lsqcurvefit(@ex08,x19,t,c);

%多次迭代，使求解更精确

A=[x1;x2;x3;x4;x5;x6;x7;x8;x9;x10;x11;x12;x13;x14;x15;x16;x17;x18;x19;x20];

B=[norm1;norm2;norm3;norm4;norm5;norm6;norm7;norm8;norm9;norm10;norm11;norm12;norm13;norm14;norm15;norm16;norm17;norm18;norm19;norm20];求得的结果如下：

ans=

4.45040.226640.861473.8909

3.99370.258344.449739.6316

3.75410.270745.810835.1695…….

3.62640.279946.787834.2331

3.62700.279946.788234.2331

3.62630.279946.788634.2330

所以，求得：k1=3.6263,k=0.2799,b=46.7886。误差的平方res=34.2330。例4BP神经网络学习算法的MATLAB实现MATLAB中BP神经网络的重要函数和基本功能newff()功能建立一个前向BP网络格式net=newff(PR，[S1S2...SN1]，{TF1TF2...TFN1}，BTF，BLF，PF)说明net为创建的新BP神经网络；PR为网络输入取向量取值范围的矩阵；[S1S2…SNl]表示网络隐含层和输出层神经元的个数；{TFlTF2…TFN1}表示网络隐含层和输出层的传输函数，默认为‘tansig’；BTF表示网络的训练函数，默认为‘trainlm’；BLF表示网络的权值学习函数，默认为‘learngdm’；PF表示性能数，默认为‘mse’。BP神经网络学习算法的MATLAB实现MATL