实验二_数值数组与遗传算法.doc

实验二 数值数组（矩阵）及其运算一、实验目的1、掌握 MATLAB 软件环境下进行数值数组（矩阵）的创建、访问的基本方法。2、掌握数值数组（矩阵）的算术运算、逻辑运算和关系运算的方法。二、实验内容1、关于“变量”（1）不需要事先定义，在遇到新的变量名时，自动创建并分配存储空间；当遇到已存在的变量时，更新其内容。（2）变量名由字母、数字和下划线构成，并且必须以字母开头。（3）能区分大小写字母。（4）对变量的赋值可采用赋值语句：变量表达式；（5）如果行末加上分号，则在屏幕上不显示结果，这在设计 M 文件时非常有用。MATLAB 中提供了一些用户不能清除的固定变量：ans，eps，pi，i，j，NaN，Inf，realmax，realmin，nargin，nargoutepsans =2.2204e-0161+eps1ans =1realmaxans =1.7977e+308realminans =2.2251e-308realmin0ans =11+realmin1ans =01+realmin=1ans =12、矩阵的创建和操作（1）矩阵的创建方法 1：输入元素列表l 矩阵行中的元素以空格（ ）或逗号（，）间隔l 矩阵行之间用分号（；）或回车（enter）间隔l 整个元素列表用方括号（ ）括起来例如: a=1 2 3;4 5 6;7 8 9a = 1 2 3 4 5 6 7 8 9l 用冒号（：）使操作简便：a=1:3;4:6;7:9a = 1 2 3 4 5 6 7 8 96方法 2：从外部数据文件中读取矩阵（见实验一中的 load 命令，略）方法 3：利用 MATLAB 内部函数产生矩阵函数功 能函数功 能compan伴随阵magic魔方矩阵zeros元素全为 0 的矩阵rand元素服从均匀分布的随机矩阵diag对角阵linspace线性空间向量ones元素全为 1 的矩阵randn元素服从正态分布的随机矩阵eye对角线上元素为 1 的矩阵（单位矩阵）logspace对数空间向量举例： b=eye(3)b = 1 0 0 0 1 0 0 0 1b = c=ones(2,5)c = 1 1 1 1 1 1 1 1 1 1d=zeros(3,2)d = 0 0 0 0 0 0 e=linspace(-3,6,10)e = -3 -2 -1 0 1 2 3 4 5 6 f=logspace(1,6,5)f = 1.0e+006 * 0.0000 0.0002 0.0032 0.0562 1.0000（2）矩阵的下标访问（存取）矩阵元素的依据全下标：行下标，列下标。如 a（i,j），几何概念清楚，引述简单。单下标（linear index）：矩阵元素在内存中是按“列”存储的“一维编号”：把矩阵的所有列按先左后右的次序首尾相接排成“一维长列”后依次编号两种下标的转换关系：A 为（mn）的二维数组（矩阵），其中某元素的“全下标”为：（r，c），则相应的“单下标”为：k=（c-1）m+rMATLAB 有两个指令可以实现以上标志方法的转换：l sub2ind %据全下标换算出单下标，如 IND = sub2ind(SIZ,I,J)l ind2sub %据单下标换算出全下标，如I,J = ind2sub(SIZ,IND) a=1 2 3; 4 5 6; 7 8 9a = 1 2 3 4 5 6 7 8 9IND = sub2ind(size(a),2,3) %返回矩阵a的2行3列元素的单下标IND = 8 aa=a(1,2)+a(2,3) %两个全下标元素值相加aa = 8 bb=a(4)+a(8) %两个单下标元素值相加，注意与上面相同。bb = 8 cc=a(4,2)? Index exceeds matrix dimensions. %访问不存在的元素时，产生错误信息a(end-2:end,2) %在下标中可以直接用 end 表示这一维的最后一个元素。 %a(end-2:end,2) 等效于a(1:3,2)ans = 2 5 8 修改元素： a(3,3)=19 %把原a(3,3)=9改成19a = 1 2 3 4 5 6 7 8 19 a(2,1:3) %显示2行1-3列ans = 4 5 6 a(2,1:3)=5 10 15 %修改成5 10 15。a = 1 2 3 5 10 15 7 8 9 扩大矩阵（增加行或列） a(4,2)=20a = 1 2 3 5 10 15 7 8 9 0 20 0（3）矩阵行列删除 a(4,:)= %删除第4行所有元素，冒号表示所有行或列a = 1 2 3 5 10 15 7 8 9 a(:,3)= %删除第3行所有元素，冒号表示所有行或列a = 1 2 5 10 7 8 a(4,2)= %删除命令应指向一个存在的元素，否则导致一个出错信息。? Indexed empty matrix assignment is not allowed. a(2:2:6)= %用单下标可删除单个或多个元素，剩余部分变成行向量。a = 1 7 10（4）矩阵连接 a=1 2;3 4; b=a a+5;a-5 zeros(size(a) %用 将小矩阵连接成大矩阵。b = 1 2 6 7 3 4 8 9 -4 -3 0 0 -2 -1 0 0 c=a;5 10c = 1 2 3 4 5 10（5）矩阵操作diag 对角矩阵和矩阵的对角化rot90 矩阵旋转 90 度reshape 矩阵元素重新排列tril 矩阵的下三角阵triu 矩阵的上三角阵fliplr 矩阵左右翻转flipud 矩阵上下翻转cat 矩阵连接repmat 复制并平铺矩阵举例：设有矩阵 a=1 2 3; 4 5 6;7 8 9a = 1 2 3 4 5 6 7 8 9 b=diag(a)b = 1 5 9% X 为矩阵时，V=diag(X,k)得到列向量 V，它取自 X 的第 k 个对角线上的元素% k=0，表示主对角线，等同于V=diag(X)。 c=diag(b) %产生对角阵。k=0，表示主对角线c = 1 0 0 0 5 0 0 0 9 c1=diag(b,1) % k0，表示在主对角线之上。c1 = 0 1 0 0 0 0 5 0 0 0 0 9 0 0 0 0 c1=diag(b,-1) % k a=1 2 3 4 %V 为 n 维向量时，diag(V,k)产生 n+abs(k)阶方阵，V 的元素处于第 k 条对角线上。a = 1 2 3 4 diag(a,2)ans = 0 0 1 0 0 0 0 0 0 2 0 0 0 0 0 0 3 0 0 0 0 0 0 4 0 0 0 0 0 0 0 0 0 0 0 0l d=rot90(a) %将矩阵 a 逆时针旋转 90 度，等同于 d=rot90(a,k), k=1 a=1 2 3 ;4 5 6;7 8 9 a = 1 2 3 4 5 6 7 8 9 d=rot90(a)d = 3 6 9 2 5 8 1 4 7l reshape(a,m,n) % reshape(a,m,n)从矩阵 a 中重新生成 mn 的矩阵，元素按“列”顺序选取。 a=1 2 3;4 5 6a = 1 2 3 4 5 6 reshape(a,3,2)ans = 1 5 4 3 2 6l a=1 2 3 ;4 5 6;7 8 9 %矩阵的下三角阵和上三角阵a = 1 2 3 4 5 6 7 8 9 f1=tril(a), f2=triu(a) f1 = 1 0 0 4 5 0 7 8 9f2 = 1 2 3 0 5 6 0 0 9l f3=fliplr(a), f4=flipud(a) %矩阵左右和上下翻转f3 = 3 2 1 6 5 4 9 8 7f4 = 7 8 9 4 5 6 1 2 3l A=1 2;3 4, B=5 6;7 8 %矩阵连接。A = 1 2 3 4/233323584?portal=homeFootprint&ref=home_footprintB = 5 6 7 8 C1=cat(1,A,B), C2=cat(2,A,B) %矩阵连接，cat（dim，A，B）将矩阵 A、B 沿着 dim 维连接起来C1 = 1 2 3 4 5 6 7 8C2 = 1 2 5 6 3 4 7 8l repmat(a,2,3) %按 mn 的格式扩张a（矩阵或标量），repmat(a,m,n) 或 repmat(a,m n)；m=n 时，repmat(a,n)a = 1 2 3 4 5 6 7 8 9 repmat(a,2,3)ans = 1 2 3 1 2 3 1 2 3 4 5 6 4 5 6 4 5 6 7 8 9 7 8 9 7 8 9 1 2 3 1 2 3 1 2 3 4 5 6 4 5 6 4 5 6 7 8 9 7 8 9 7 8 9 repmat(8,3)ans = 8 8 8 8 8 8 8 8 83、矩阵运算（1）算术运算1）算术运算符： * / 2）算术运算和对应的 M 函数功 能算术表达式M 函数功 能算术表达M 函数加 法A+Bplus(A,B)矩阵左除A Bmldivide(A,B)减 法A-Bminus(A,B)数组左除A Bldivide(A,B)矩阵乘法A*Bmtimes(A,B)矩阵求幂ABmpower(A,B)数组乘法A*Btimes(A,B)数组求幂ABpower(A,B)矩阵右除A/Bmrdivide(A,B)矩阵共轭转置A ctranspose(A)数组右除A/Brdivide(A,B)非共轭转置A transpose(A)3）基本数学函数： 函数名 数 学 计 算 功 能 函 数 名 数 学 计 算 功 能abs（x）实数的绝对值或复数的幅值floor（x）对 x 朝-方向取整acos（x）反余弦 arcsin xgcd（m，n）求正整数 m 和 n 的最大公约数acosh（x）反双曲余弦 arccosh ximag（x）求复数 x 的虚部angle（x）在四象限内求复数 x 的相角lcm（m，n）求正整数 m 和 n 的最小公倍数asin（x）反正弦 arcsin xlog（x）自然对数（以e 为底数）asinh（x）反双曲正弦 arcsinh xlog10（x）常用对数（以 10 为底数）atan（x）反正切 arctan xreal（x）求复数 x 的实部atan2（x,y）在四象限内求反正切rem（m，n）求正整数 m 和 n 的 m/n 之余数atanh（x）反双曲正切 arctanh xround（x）对 x 四舍五入到最接近的整数ceil（x）对 x 朝+方向取整sign（x）符号函数：求出 x 的符号conj（x）求复数 x 的共轭复数sin（x）正弦 sin xcos（x）余弦 cos xsinh（x）反双曲正弦 sinh xcosh（x）双曲余弦 cosh xsqrt（x）求实数 x 的平方根exp（x）指数函数 extan（x）正切 tan xfix（x）对 x 朝原点方向取整tanh（x）双曲正切 tanh x举例：两种不同转置的区别：clear; A=zeros(2,3);A(:)=1:6; %全元素赋值法A=A*(1+i) %运用标量与数组乘产生复数矩阵A_A=A. %数组转置，即非共轭转置A_M=A %矩阵转置，即共轭转置A = %运用标量与数组乘产生复数矩阵。 1.0000 + 1.0000i 3.0000 + 3.0000i 5.0000 + 5.0000i 2.0000 + 2.0000i 4.0000 + 4.0000i 6.0000 + 6.0000iA_A = %数组转置，即非共轭转置。 1.0000 + 1.0000i 2.0000 + 2.0000i 3.0000 + 3.0000i 4.0000 + 4.0000i 5.0000 + 5.0000i 6.0000 + 6.0000iA_M = %矩阵转置，即共轭转置1.0000 - 1.0000i 2.0000 - 2.0000i 3.0000 - 3.0000i 4.0000 - 4.0000i 5.0000 - 5.0000i 6.0000 - 6.0000i求矩阵的逆矩阵（4 种方法）：B=1 2;3 4;B1=inv(B)B2=B(-1)B3=eye(2)/BB4=Beye(2)B1 = -2.0000 1.0000 1.5000 -0.5000B2 = -2.0000 1.0000 1.5000 -0.5000B3 = -2.0000 1.0000 1.5000 -0.5000B4 = -2.0000 1.00001.5000 -0.5000关于矩阵求幂：A=1 2;3 4;B=2 1;3 2;A.Bans = 1 227 16exp(A) %求eA ，其他函数如 sin，tan， log，sqrt 等都是对矩阵的每个元素进行相应运算。ans = 2.7183 7.3891 20.0855 54.5982mpower(A,B) %函数操作数 A、B 中至少要有一个标量；否则出错。? Error using = mpowerAt least one operand must be scalar. B=5;mpower(A,B)ans = 1069 1558 2337 3406 A(-2)ans = 5.5000 -2.5000 -3.7500 1.7500%Xp，标量 p 为整数时，通过连乘实现；当 p为负整数时，先对矩阵X求逆。%当 p 为其他值时，计算要涉及到特征值和特征向量，这样当V，D=eig（X）Xp=V*D.p/V（2）关系和逻辑运算1）关系操作符（6 种）：、=、all(a) %a 为向量时，all（a）在所有元素非零时，得到逻辑真：“1”%a 为矩阵时，all（a）将 a 的列当作向量，得到 0、1 构成的行向量。ans =0 1 all(all(a) %对矩阵 a 采用两次 all 函数可得到标量结果。ans = 0 all(ones(2,5),2) %all（a,dim）可测试沿着 dim 维上的逻辑条件。ans = 1 1any(a) %any：测试“任意元素”是否为非零值，也有 any（a）、any（a，dim）等形式。ans =1 1a=zeros(5,10) %find：查找非零元素的值及其下标，k=find（a）：找非零元素的下标，若a全零，则 k= a=zeros(5,10)a = 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 a(3,7)=0.5a = 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0.500 0 0 00 0 0 0 0 0 0 0 0 0 k=find(a) % 找非零元素的单下标。k =33 i,j=find(a) %找出矩阵 a 中非零元素的行列下标。i = 3j = 7 i,j,v=find(a) %i,j,v=find（X）：找出矩阵 X 中非零元素的行列下标和值。i = 3j = 7v = 0.5000% %后两者的功能常用于“稀疏矩阵”中。l if exist(sg.dat) %exist：可在装入数据之前对数据文件作检测（存在否？）load sg.mat %这在数据量较大时是很有用的else %sg.dat 存在时直接将数据读入，否则将 sg 初始化为 0 矩阵sg=zeros(30,2)endl 利用 is*（如：isreal、ischar、islogical 等）这一组函数可以对矩阵进行各种检测，下面以 isnan(用于检测非数值，同 isinf)举例：m=mean(a) %求均值；矩阵中含 NaN 时，基于它的函数值也含有 NaN；因此，在数据处理之前，一般应对数据进行分析，删除含有 NaN 的测量样本后再进行处理。 a=1 2 2;2 3 0.4; 3 5 -0.9;4 NaN 0.9;5 -2.0 -0.8m=mean(a)i,j=find(b) if length(i) %当 i 不为空时a(i,:)= enddisp(The data after processing is a=),disp(a)a = 1.0000 2.0000 2.0000 2.0000 3.0000 0.4000 3.0000 5.0000 -0.9000 4.0000 NaN 0.9000 5.0000 -2.0000 -0.8000m = 3.0000 NaN 0.3200i = 4j = 1a = 1.0000 2.0000 2.0000 2.0000 3.0000 0.4000 3.0000 5.0000 -0.9000 5.0000 -2.0000 -0.8000The data after processing is a= 1.0000 2.0000 2.0000 2.0000 3.0000 0.4000 3.0000 5.0000 -0.90005.0000 -2.0000 -0.80004. 遗传算法的Matlab程序%function test(lchrom,popsize,Pc,Pm,gen)%function test(xy,lchrom,popsize,Pc,Pm,gen) %初始化数据lchrom=10 %基因个数=城市个数popsize=100 %种群规模Pc=0.91 %交叉概率Pm=0.08 %变异概率gen=1000 %最大迭代次数tic %计时开始%xy=rand(lchrom,2)*100 %随机给出一组初始点for i=1:lchrom plot(xy(i,1),xy(i,2),r*); hold onend%55shortpath= ; % 最短路径设一个最大值shortpathlength=realmax;pop=inigroup(lchrom,popsize) %给出一组popsize大小的初始化排序方案record=zeros(1,gen)for w=1:gen %for 1 %最大迭代次数% for R=1:popsize %for 2 %种群规模，一一计算各染色体 if rand(1)max(t3(2),t4(2) %if 2 %保留交叉后四个染色体中较好的染色体，适应值越大越好，因为适应值=1/距离 pop(i,:)=oldp1; pop(j,:)=oldp2; if t1(2)t2(2) %if 3 goodp=oldp1; t=t1; else goodp=oldp2; t=t2; end %end if 3 else pop(i,:)=newp1; pop(j,:)=newp2; if t3(2)t4(2) %if 4 goodp=newp1; t=t3; else goodp =newp2; t=t4; end %end if 4 end %end if 2 if t(1)shortpathlength %if 5 %保留当前最佳适应值 shortpath=goodp; shortpathlength=t(1); end %end if 5 end %end if 1 %完成一次染色体交叉， end %end for 1 随机交叉最多popsize次% for R=1:popsize %for 3 %实现变异 if rand(1)Pm %if 6 %满足变异概率才进行变异 p=pop(R,:); %选择产生变异的染色体 a=floor(lchrom*rand(1)+1; %产生1-10随机数 b=floor(lchrom*rand(1)+1; minmutation=max(a,b); % 例如，minmutation=10，maxmutation=7 maxmutation=min(a,b); %pop3,:=6 10 2 9 7 1 4 8 5 3第3个染色体变异 p(minmutation maxmutation)=p(maxmutation minmutation); %第7和10个基因交换变异，这样可避免产生非法变异 t=pathlenfit(p,distmatrix); %计算变异后的适应值 if t(1)shortpathlength %if 7 如果变异后最后，保留之。 shortpath=p; shortpathlength=t(1); end %end if 7 end %end if 6 end %end for 3 %最多实现popsize次变异 % for R=1:popsize %for 4 % 逆序操作，增加方案的多样性 if rand(1)2*Pm %if 8 a=floor(lchrom-1)*rand(1)+1; b=floor(lchrom-1)*rand(1)+1; mininverse=min(a,b); % 例如，mininverse=4；maxinverse=9 maxinverse=max(a,b); order=1:lchrom; % order=1 2 3 4 5 6 7 8 9 10 order(mininverse:maxinverse)=order(maxinverse:-1:mininverse); %从4到9逆序排列 order=1 2 3 9 8 7 6 5 4 10 pold=pop(R,:); %R=18，该染色体被选中，设当前顺序是6 10 1 4 2 8 5 3 9 7 told=pathlenfit(pold,distmatrix); %计算旧排列适应值 pnew=pold(order); %按照order给定顺序重排染色体6 10 1 9 3 5 8 2 4 7 tnew=pathlenfit(pnew,distmatrix); %计算新排列适应值 if tnew(1)told(1) %if 9 保留当前最好适应值 pop(R,:)=pnew; if tnew(1)shortpathlength %if 10 shortpath=pnew; shortpathlength=tnew(1); end %end if 10 end %end if 9 end %end if 8 end %end for 4%record(w)=shortpathlength; %记录当前迭代次数下最佳适应值end %end for 1 %达到循环迭代次数draw(shortpath,xy) %画出最终路径图。shortpath %打印最短路径。shortpathlengthtoc %打印计算时间function pop=inigroup(lchrom,popsize)% lchrom-城市数% popsize-初始种群数% pop-初始种群 行数为popsize ,列数为lchrom n=1; pop=zeros(popsize,lchrom); while n=popsize pop(n,:)=randperm(lchrom); n=n+1; endfunction m=crossover(oldp1,oldp2)% oldp1,oldp2 杂交，产生m=newp1;newp2 crossj1=floor(length(oldp1)-1)*rand)+1; crossj2=floo