02流水线车间生产调度的遗传算法MATLAB源代码
- 格式:docx
- 大小:19.42 KB
- 文档页数:6
function [Zp,Y1p,Y2p,Y3p,Xp,LC1,LC2]=JSPGA(M,N,Pm,T,P)%greemsim原创--------------------------------------------------------------------------% JSPGA.m% 输入参数列表% M 遗传进化迭代次数% N 种群规模(取偶数)% Pm 变异概率% T m×n的矩阵,存储m个工件n个工序的加工时间% P 1×n的向量,n个工序中,每一个工序所具有的机床数目% 输出参数列表% Zp 最优的Makespan值% Y1p 最优方案中,各工件各工序的开始时刻,可根据它绘出甘特图% Y2p 最优方案中,各工件各工序的结束时刻,可根据它绘出甘特图% Y3p 最优方案中,各工件各工序使用的机器编号% Xp 最优决策变量的值,决策变量是一个实数编码的m×n矩阵% LC1 收敛曲线1,各代最优个体适应值的记录% LC2 收敛曲线2,各代群体平均适应值的记录% 最后,程序还将绘出三副图片:两条收敛曲线图和甘特图(各工件的调度时序图)%第一步:变量初始化[m,n]=size(T);%m是总工件数,n是总工序数Xp=zeros(m,n);%最优决策变量LC1=zeros(1,M);%收敛曲线1LC2=zeros(1,N);%收敛曲线2%第二步:随机产生初始种群farm=cell(1,N);%采用细胞结构存储种群for k=1:NX=zeros(m,n);for j=1:nfor i=1:mX(i,j)=1+(P(j)-eps)*rand;endendfarm{k}=X;endcounter=0;%设置迭代计数器while counter<M%停止条件为达到最大迭代次数? %第三步:交叉newfarm=cell(1,N);%交叉产生的新种群存在其中Ser=randperm(N);for i=1:2:(N-1)A=farm{Ser(i)}; %父代个体B=farm{Ser(i+1)};Manner=unidrnd(2);%随机选择交叉方式if Manner==1cp=unidrnd(m-1);%随机选择交叉点%双亲双子单点交叉a=[A(1:cp,:);B((cp+1):m,:)];%子代个体b=[B(1:cp,:);A((cp+1):m,:)];elsecp=unidrnd(n-1);%随机选择交叉点a=[A(:,1:cp),B(:,(cp+1):n)];b=[B(:,1:cp),A(:,(cp+1):n)];endnewfarm{i}=a;%交叉后的子代存入newfarmnewfarm{i+1}=b;end%新旧种群合并FARM=[farm,newfarm];?%第四步:选择复制FITNESS=zeros(1,2*N);fitness=zeros(1,N);plotif=0;for i=1:(2*N)X=FARM{i};Z=COST(X,T,P,plotif);%调用计算费用的子函数FITNESS(i)=Z;end%选择复制采取两两随机配对竞争的方式,具有保留最优个体的能力Ser=randperm(2*N);for i=1:Nf1=FITNESS(Ser(2*i-1));?f2=FITNESS(Ser(2*i));if f1<=f2farm{i}=FARM{Ser(2*i-1)};fitness(i)=FITNESS(Ser(2*i-1));elsefarm{i}=FARM{Ser(2*i)};fitness(i)=FITNESS(Ser(2*i));endend%记录最佳个体和收敛曲线minfitness=min(fitness)meanfitness=mean(fitness)LC1(counter+1)=minfitness;%收敛曲线1,各代最优个体适应值的记录LC2(counter+1)=meanfitness;%收敛曲线2,各代群体平均适应值的记录pos=find(fitness==minfitness);Xp=farm{pos(1)};?%第五步:变异for i=1:Nif Pm>rand;%变异概率为PmX=farm{i};I=unidrnd(m);J=unidrnd(n);X(I,J)=1+(P(J)-eps)*rand;farm{i}=X;endendfarm{pos(1)}=Xp;?counter=counter+1end%输出结果并绘图figure(1);plotif=1;X=Xp;[Zp,Y1p,Y2p,Y3p]=COST(X,T,P,plotif);figure(2);plot(LC1);figure(3);plot(LC2);function [Zp,Y1p,Y2p,Y3p]=COST(X,T,P,plotif)% JSPGA的内联子函数,用于求调度方案的Makespan值% 输入参数列表% X 调度方案的编码矩阵,是一个实数编码的m×n矩阵% T m×n的矩阵,存储m个工件n个工序的加工时间% P 1×n的向量,n个工序中,每一个工序所具有的机床数目% plotif 是否绘甘特图的控制参数% 输出参数列表% Zp 最优的Makespan值% Y1p 最优方案中,各工件各工序的开始时刻% Y2p 最优方案中,各工件各工序的结束时刻% Y3p 最优方案中,各工件各工序使用的机器编号%第一步:变量初始化[m,n]=size(X);Y1p=zeros(m,n);Y2p=zeros(m,n);Y3p=zeros(m,n);%第二步:计算第一道工序的安排Q1=zeros(m,1);Q2=zeros(m,1);R=X(:,1);%取出第一道工序Q3=floor(R);%向下取整即得到各工件在第一道工序使用的机器的编号%下面计算各工件第一道工序的开始时刻和结束时刻for i=1:P(1)%取出机器编号pos=find(Q3==i);%取出使用编号为i的机器为其加工的工件的编号lenpos=length(pos);if lenpos>=1Q1(pos(1))=0;Q2(pos(1))=T(pos(1),1);if lenpos>=2for j=2:lenposQ1(pos(j))=Q2(pos(j-1));Q2(pos(j))=Q2(pos(j-1))+T(pos(j),1);endendendendY1p(:,1)=Q1;Y2p(:,1)=Q2;Y3p(:,1)=Q3;%第三步:计算剩余工序的安排for k=2:nR=X(:,k);%取出第k道工序Q3=floor(R);%向下取整即得到各工件在第k道工序使用的机器的编号%下面计算各工件第k道工序的开始时刻和结束时刻for i=1:P(k)%取出机器编号pos=find(Q3==i);%取出使用编号为i的机器为其加工的工件的编号lenpos=length(pos);if lenpos>=1EndTime=Y2p(pos,k-1);%取出这些机器在上一个工序中的结束时刻POS=zeros(1,lenpos);%上一个工序完成时间由早到晚的排序for jj=1:lenposMinEndTime=min(EndTime);ppp=find(EndTime==MinEndTime);POS(jj)=ppp(1);EndTime(ppp(1))=Inf;end?%根据上一个工序完成时刻的早晚,计算各工件第k道工序的开始时刻和结束时刻Q1(pos(POS(1)))=Y2p(pos(POS(1)),k-1);Q2(pos(POS(1)))=Q1(pos(POS(1)))+T(pos(POS(1)),k);%前一个工件的结束时刻if lenpos>=2for j=2:lenposQ1(pos(POS(j)))=Y2p(pos(POS(j)),k-1);%预定的开始时刻为上一个工序的结束时刻if Q1(pos(POS(j)))<Q2(pos(POS(j-1)))%如果比前面的工件的结束时刻还早Q1(pos(POS(j)))=Q2(pos(POS(j-1)));endendendendendY1p(:,k)=Q1;Y2p(:,k)=Q2;Y3p(:,k)=Q3;end%第四步:计算最优的Makespan值Y2m=Y2p(:,n);Zp=max(Y2m);%第五步:绘甘特图if plotiffor i=1:mfor j=1:nmPoint1=Y1p(i,j);mPoint2=Y2p(i,j);mText=m+1-i;PlotRec(mPoint1,mPoint2,mText);Word=num2str(Y3p(i,j));text(0.5*mPoint1+0.5*mPoint2,mText-0.5,Word); hold onx1=mPoint1;y1=mText-1;x2=mPoint2;y2=mText-1;x3=mPoint2;y3=mText;x4=mPoint1;y4=mText;fill([x1,x2,x3,x4],[y1,y2,y3,y4],'r');fill([x1,x2,x3,x4],[y1,y2,y3,y4],[1,0.5,1]);text(0.5*mPoint1+0.5*mPoint2,mText-0.5,Word); endendendfunction PlotRec(mPoint1,mPoint2,mText)% 此函数画出小矩形% 输入:?% mPoint1 输入点1,较小,横坐标% mPoint2 输入点2,较大,横坐标% mText 输入的文本,序号,纵坐标vPoint = zeros(4,2) ;vPoint(1,:) = [mPoint1,mText-1];vPoint(2,:) = [mPoint2,mText-1];vPoint(3,:) = [mPoint1,mText];vPoint(4,:) = [mPoint2,mText];plot([vPoint(1,1),vPoint(2,1)],[vPoint(1,2),vPoint(2,2)]); hold on ;plot([vPoint(1,1),vPoint(3,1)],[vPoint(1,2),vPoint(3,2)]); plot([vPoint(2,1),vPoint(4,1)],[vPoint(2,2),vPoint(4,2)]); plot([vPoint(3,1),vPoint(4,1)],[vPoint(3,2),vPoint(4,2)]);。
% f(x)=11*sin(6x)+7*cos(5x),0<=x<=2*pi;%%初始化参数L=16;%编码为16位二进制数N=32;%初始种群规模M=48;%M个中间体,运用算子选择出M/2对母体,进行交叉;对M个中间体进行变异T=100;%进化代数Pc=0.8;%交叉概率Pm=0.03;%%变异概率%%将十进制编码成16位的二进制,再将16位的二进制转成格雷码for i=1:1:Nx1(1,i)= rand()*2*pi;x2(1,i)= uint16(x1(1,i)/(2*pi)*65535);grayCode(i,:)=num2gray(x2(1,i),L);end%% 开始遗传算子操作for t=1:1:Ty1=11*sin(6*x1)+7*cos(5*x1);for i=1:1:M/2[a,b]=min(y1);%找到y1中的最小值a,及其对应的编号bgrayCodeNew(i,:)=grayCode(b,:);%将找到的最小数放到grayCodeNew中grayCodeNew(i+M/2,:)=grayCode(b,:);%与上面相同就可以有M/2对格雷码可以作为母体y1(1,b)=inf;%用来排除已找到的最小值endfor i=1:1:M/2p=unidrnd(L);%生成一个大于零小于L的数,用于下面进行交叉的位置if rand()<Pc % Pc是交叉概率%将选定的染色体的点后的基因进行交换for j=p:1:Ltemp= grayCodeNew(i,j);grayCodeNew(i,j)=grayCodeNew(M-i+1,j);grayCodeNew(M-i+1,j)=temp;endendend%%将全部染色体按概率进行变异for i=1:1:Mfor j=1:1:Lif rand()<Pm %Pm为染色体变异的概率grayCodeNew(i,j)=dec2bin(1-bin2dec(grayCodeNew(i,j)));endendend%%第一代结束后生成的较优的染色体,得以保存然后进行下一代操作for i=1:1:Mx4(1,i)=gray2num(grayCodeNew(i,:));endx3=double(x4)*2*pi/65535;y3=11*sin(6*x3)+7*cos(5*x3);for i=1:1:N[a,b]=min(y3);x1(1,i)=x3(1,b);grayCode(i,:)=grayCodeNew(b,:);y3(1,b)=inf;endendx1y1=11*sin(6*x1)+7*cos(5*x1)。
$标准遗传算法%优化函数为f=-(x-l厂2+4,其中,0<=x<=3$编码长度为10位,编码精度为0. 0029嗚种群规模设为40,遗传算J'•分别为比例选择,单点交叉和单点变异。
交叉槪率0・7,变异概率0・1 %最大进化代数为200代,保优操作。
main・ minitial;global G;for G=l:200crossover;mutation;adapting;keepbest; selection;endresult;嗚初始化函数,甌机形成规模为40初始种群initial.mpop (40, 10)=0;best_individual (10)=0; %最优个体adaptive (200) =0; %种群平均适应值for i=l:40for j=l:10if rand>0・ 5pop(i, j)=l;elsepop(i, j)=0;endendend% popclear i;clear j;务交叉操作,概率为0.7,单点交叉crossover・ mfor i=l:2:39cross_P=rand; $随机产生•个数,以比较交叉概率if cross_P<0. 9 $交叉概率为0・9 cross_pos=round(10*rand);玄交叉位置为0为,若位置为0或9,则不辿〔if or (cross_pos=0, cross_pos==9)continue;for j=cross_pos:10 temp二pop (i, j); pop(i, j)=pop(i+l, j); pop(i+l, j)=temp;endendendclear i;clear j;clear temp;clear cross_P;clear cross_pos;賀变异操作,单点变异,变异概率为0.1 mutation・ mfor i=l:40if rand<0. 1 $通过变异槪率M_pos=round(l0*rand);if Mj)os"=0 %若变异位为0则无总义pop(i, M_pos)=l-pop(i, M_pos);endendendclear M_pos;clear i;%计算适应值adapting・ mfor i=l:40adapt(i)=0;endfor i=l:40for j=l:10if pop(i, j)==l adapt (i)=adapt (i)+2* (10-j);endendadapt (i)=adapt (i)*0. 0029; adapt (i)=-(adapt (i)-l). *2+4;endglobal adapt_best;global best_pos;adapt_best=0; %最佳个体best_pos=0; %最佳个体在种群「11的位置% adapt_ave=O;for i=l:40adaptive(G) =adapt__avc(G)-t-adapt(i):i f adapt__best<adapt(i) adapt_bcst=adapt(i); best_pos=i;endendadaptive (G) =adapt_ave (G)/ 40;clear i:clear j;制呆优操作keepbest・ mfor i=l:10best indi v idual(i)=pop(best pos, i);end% The select oprator functionselection・ mada_sum=0;ada^temp^O;r=0;i=0;j=0;for i=l:40ada_sum=ada_sum+adapt(i):endfor i=l:39 %选抒39次,M.n; -个个体留给历代最优解r rand*ada_sum; %随机产生_个数 ada.temp=0: %初始化累加值为0 j=0;whi le(ada_temp<r)j=j+l;ada_temp=ada_temp+adapt(j):end%退出循环时的j值即为被选择的个体序号% if j>40% j=40;% endfor k=l:10new_pop(i, k)=pop(j, k);end %最优解复制for i=l:10new_pop(40, i)=best_individual (i); end玄将选择产生的新群体复制给pop种群for i=l:40for j=l:10pop(i, j)=new_pop(i, j);endendclear i;clear j;clear k;clear r;clear ada_temp;賀结果统计函数result・ mplot(adapt_ave);best=0;for j=l:10if best_individual(j)==1 best=best+2*(10-j);endendbest=best*O. 0029;'最优个体为'best'最优解为'best=-(best-1). "2+4;best。
matlab鸟群算法求解车间调度问题详解及实现源码⽬录⼀、车间调度简介1 车间调度定义2 传统作业车间调度3 柔性作业车间调度⼆、蝴蝶优化算法(MBO)简介1 介绍2 ⾹味3 具体算法三、部分源代码五、matlab版本及参考⽂献⼀、车间调度简介1 车间调度定义车间调度是指根据产品制造的合理需求分配加⼯车间顺序,从⽽达到合理利⽤产品制造资源、提⾼企业经济效益的⽬的。
车间调度问题从数学上可以描述为有n个待加⼯的零件要在m台机器上加⼯。
问题需要满⾜的条件包括每个零件的各道⼯序使⽤每台机器不多于1次,每个零件都按照⼀定的顺序进⾏加⼯。
2 传统作业车间调度传统作业车间带调度实例有若⼲⼯件,每个⼯件有若⼲⼯序,有多个加⼯机器,但是每道⼯序只能在⼀台机器上加⼯。
对应到上⾯表格中的实例就是,两个⼯件,⼯件J1有三道⼯序,⼯序Q11只能在M3上加⼯,加⼯时间是5⼩时。
约束是对于⼀个⼯件来说,⼯序的相对顺序不能变。
O11->O12->O13。
每时刻,每个⼯件只能在⼀台机器上加⼯;每个机器上只能有⼀个⼯件。
调度的任务则是安排出⼯序的加⼯顺序,加⼯顺序确定了,因为每道⼯序只有⼀台机器可⽤,加⼯的机器也就确定了。
调度的⽬的是总的完⼯时间最短(也可以是其他⽬标)。
举个例⼦,⽐如确定了O21->O22->O11->O23->O12->O13的加⼯顺序之后,我们就可以根据加⼯机器的约束,计算出总的加⼯时间。
M2加⼯O21消耗6⼩时,⼯件J2当前加⼯时间6⼩时。
M1加⼯O22消耗9⼩时,⼯件J2当前加⼯时间6+9=15⼩时。
M3加⼯O11消耗5⼩时,⼯件J1当前加⼯时间5⼩时。
M4加⼯O23消耗7⼩时,⼯件J2加⼯时间15+7=22⼩时。
M1加⼯O12消耗11⼩时,但是要等M1加⼯完O22之后才开始加⼯O12,所以⼯件J1的当前加⼯时间为max(5,9)+11=20⼩时。
M5加⼯O13消耗8⼩时,⼯件J2加⼯时间20+8=28⼩时。
01.02.% P MX me ans G oldbe rg'sParti allyMappe d Cro Ssove r)03.% P roced ure :PMX04.%Step1. Sel ect t wo po sitio ns al ong t he st ringunifo rmlyat ra ndom.05.% T he su bstri ngs d efine d bythe t wo po sitio ns ar e cal led t he ma pping sect ions. 06.% N ote:(You c an wr ite t o sel ect a star t poi nt an d a l ength)07.% S tep2. Exch angetwo s ubstr ingsbetwe en pa rents to p roduc e pro to-ch ildre n.08.% S tep3. Dete rmine themappi ng re latio nship betw een t wo ma pping sect ion.09.% Step4. Le galiz e off sprin g wit h the mapp ing r elati onshi p.10.fun ction [new Va,ne wVb]=PMX1(Va,Vb)11.fpri ntf('origi nal V a and Vb a re:\n')12.%Va= [ 1 6 10 3 9 4 52 7 8 ];13.%V b= [2 9 1 4 10 5 6 8 37 ];14.%Va=1:915.%Vb=[5 46 9 2 1 78 3]16.%--------------------------------------------------------------------------------17.%Step1. Se lecttwo p ositi ons a longthe s tring unif ormly at r andom.18.star tXorP oint=mod(c eil(r and(1)*10),leng th(Va) );19.i f sta rtXor Point==020. star tXorP oint=start XorPo int+1;21.end22.xor Lengt h=mod(floo r(ran d(1)*10),l ength(Va));23.endX orPoi nt=st artXo rPoin t+xor Lengt h;24.whi le(en dXorP oint>lengt h(Vb) )25. xorLe ngth=mod(f loor(rand(1)*10),len gth(V a));26. end XorPo int=s tartX orPoi nt+xo rLeng th;27.en d 28.f print f('\n The(star tXorP oint,endXo rPoin t)=(%d,%d)\n',s tartX orPoi nt,en dXorP oint) 29.%star tXorP oint=330.%end XorPo int=631.%--------------------------------------------------------------------------------32.% St ep2.Excha nge t wo su bstri ngs b etwee n par entsto pr oduce prot o-chi ldren.33.temp1=Va(start XorPo int:e ndXor Point);34.tem p2=Vb(star tXorP oint:endXo rPoin t);35.Va(star tXorP oint:endXo rPoin t)=te mp2;36.V b(sta rtXor Point:endX orPoi nt)=t emp1;37.clear temp1;38.cle ar te mp2;39.f print f('Th e exc hange d Vaand V b are:\n')40.Va41.Vb42.%--------------------------------------------------------------------------------43.% Ste p3. D eterm ine t he ma pping rela tions hip b etwee n two mapp ing s ectio n.44.tem p1=Va(star tXorP oint:endXo rPoin t);45.te mp2=V b(sta rtXor Point:endX orPoi nt);46.f or ix=1:le ngth(temp1)47. rawMa pRela tion(ix,1:2)=[V a(sta rtXor Point+ix-1),Vb(start XorPo int+i x-1)];48.end49.rawM apRel ation50.%rawM apRel ation=[6 3;9 4;2 5;1 6;37]51.row Index=1;52.co lInde x=1;53.w hile( rowI ndex<=size(rawM apRel ation,1) )54. wh ile(colIn dex<=size(rawMa pRela tion,2) )55.rawMa pRela tion(rowIn dex,c olInd ex )56. [i,j]=fin d(raw MapRe latio n==ra wMapR elati on(ro wInde x,col Index ) ) ;57. if(leng th(i)>1)58. if( j(1)<j(2) )59. t empRe sult=[rawM apRel ation(i(2),:),rawMa pRela tion(i(1),:)];60. k=161. whil e k<l ength(temp Resul t)62. if temp Resul t(1,k)==te mpRes ult(1,k+1)63. temp Resul t(k:k+1)=[];64. en d 65.k=k+1;66.end67. tem pResu lt 68. raw MapRe latio n(i,:)=[];69. ra wMapR elati on(si ze(ra wMapR elati on,1)+1,1:2)=te mpRes ult;70. 71.e lse 72. tem pResu lt=[r awMap Relat ion(i(1),:), ra wMapR elati on(i(2),:)]; 73. k=174.w hilek<len gth(t empRe sult)75. if t empRe sult(1,k)==temp Resul t(1,k+1)76. t empRe sult(k:k+1)=[];77. end78. k=k+1;79. en d80. tempR esult81. raw MapRe latio n(i,:)=[];82. ra wMapR elati on(si ze(ra wMapR elati on,1)+1,1:2)=te mpRes ult;83. 84. end85. end86. if(lengt h(i)==1 &lengt h(j)==1) 87. col Index=colI ndex+1;88. els e89. r owInd ex=190. col Index=1;91. en d 92. end93. row Index=rowI ndex+1;94. end 95. col Index=1;%R eset96.97.rawM apRel ation98.tMap=[rawM apRel ation;flip lr(ra wMapR elati on)]99.M ap=tM ap'100.f print f('\n The(star tXorP oint,endXo rPoin t)=(%d,%d)\n',s tartX orPoi nt,en dXorP oint) 101.Va102.V b103.%--------------------------------------------------------------------------------104.% Step4. Le galiz e off sprin g wit h the mapp ing r elati onshi p.105.if star tXorP oint~=1106. fori=1:s tartX orPoi nt-1107. [r,c]=fi nd(Ma p(1,:)==Va(1,i)) ; 108. if~isem pty(r) & ~isemp ty(c)109. Va(1,i)=Ma p(r+1,c);110. end111. [r1,c1]=fin d(Map(1,:)==Vb(1,i)); 112. if ~isemp ty(r1) & ~isemp ty(c1)113. Vb(1,i)=M ap(r1+1,c1);114. e nd 115. e nd116.en d117.ifendXo rPoin t~=le ngth(Va)118. for i=en dXorP oint+1:len gth(V a)119. [r,c]=find(Map(1,:)==Va(1,i));120. i f ~is empty(r) & ~ise mpty(c)121. Va(1,i)=Map(r+1,c);122. en d 123. [r1,c1]=f ind(M ap(1,:)==V b(1,i)) ;124. if ~ise mpty(r1) & ~ise mpty(c1) 125. Vb(1,i)=Map(r1+1,c1);126. end127. end128.en d 129.fprin tf('T he fi nal V a and Vb a re:\n') 130.ne wVa=V a131.new Vb=Vb132.。
遗传算法优化相关MATLAB算法实现遗传算法(Genetic Algorithm,GA)是一种基于生物进化过程的优化算法,能够在空间中找到最优解或接近最优解。
它模拟了自然选择、交叉和变异等进化操作,通过不断迭代的方式寻找最佳的解。
遗传算法的主要步骤包括:初始化种群、评估适应度、选择、交叉、变异和更新种群等。
在MATLAB中,可以使用遗传算法工具箱(Genetic Algorithm & Direct Search Toolbox)来实现遗传算法的优化。
下面以实现一个简单的函数优化为例进行说明。
假设我们要优化以下函数:```f(x)=x^2-2x+1```首先,我们需要定义适应度函数,即上述函数f(x)。
在MATLAB中,可以使用如下代码定义适应度函数:```MATLABfunction fitness = myFitness(x)fitness = x^2 - 2*x + 1;end```接下来,我们需要自定义遗传算法的参数,包括种群大小、迭代次数、交叉概率和变异概率等。
在MATLAB中,可以使用如下代码定义参数:```MATLABpopulationSize = 100; % 种群大小maxGenerations = 100; % 迭代次数crossoverProbability = 0.8; % 交叉概率mutationProbability = 0.02; % 变异概率```然后,我们需要定义遗传算法的上下界范围。
在本例中,x的范围为[0,10]。
我们可以使用如下代码定义范围:```MATLABlowerBound = 0; % 下界upperBound = 10; % 上界```接下来,我们可以使用遗传算法工具箱中的`ga`函数进行遗传算法的优化。
如下所示:```MATLAB```最后,我们可以得到最优解x和最优值fval。
在本例中,我们得到的结果应该接近1以上只是一个简单的例子,实际应用中可能需要根据具体问题进行参数的设定和函数的定义。
%用遗传算法求函数f(x1,x2)=100(x1^2-x2)^2+(1-x1)^2的最大值,-2.048<=x1,x2<=2.048 Umax=2.048;Umin=-2.048;N=80;%种群规模T=100;%迭代到第几代Pc=0.8;%交叉概率Pm=0.05;%变异概率L=10;%编码长度bval=round(rand(N,2*L));%随机产生初始种群的80个个体obj=zeros(1,N);%种群每个个体对应的函数f(x1,x2)的值for ii=1:Tfor i=1:Ny1=0;y2=0;for j=1:Ly1=y1+bval(i,L-j+1)*2^(j-1);endfor j=1:Ly2=y2+bval(i,2*L-j+1)*2^(j-1);end%把种群中个体解码为对应的x1和x2x1=(Umax-Umin)*y1/(2^L-1)+Umin;%把编码左边L位解码为x1x2=(Umax-Umin)*y2/(2^L-1)+Umin;%把编码右边L位解码为x2obj(i)=100*(x1*x1-x2).^2+(1-x1).^2;%计算x1和x2对应的函数f(x1,x2)的值,即每个个体的适应度xx(i,:)=[x1,x2];endfunc=obj;%适应度函数%fsum=sum(func);p=func/sum(func);%求每个个体的选择概率q=cumsum(p);%求累加概率[fmax,indmax]=max(func);%求最佳个体bestv=-inf;if fmax>=bestvbestv=fmax;bvalxx=bval(indmax,:);optxx=xx(indmax,:);endfor i=1:N-1r=rand;tmp=find(r<=q);newbval(i,:)=bval(tmp(1),:);endnewbval(N,:)=bvalxx;%最优个体保留bval=newbval;for i=1:2:(N-1)cc=rand;%单点交叉if cc<Pcpoint=ceil(rand*(2*L-1));%取一个1到2L-1的整数,作为交叉点ch=bval(i,:);bval(i,point+1:2*L)=bval(i+1,point+1:2*L);bval(i+1,point+1:2*L)=ch(1,point+1:2*L);end%变异cm=rand;if cm<Pmpoint=ceil(rand*(2*L-1));%取一个1到2L-1的整数,作为变异点bval(i,point)=~bval(i,point);endendendbestv,optxx%绘图,观察函数图象,确定最大值点的位置[x,y]=meshgrid(-2.05:0.01:2.05);z=100*(x.*x-y).^2+(1-x).^2;mesh(x,y,z);x1=2.048,y1=2.048,z1=100*(x1.*x1-y1).^2+(1-x1).^2x2=2.048,y2=-2.048z2=100*(x2.*x2-y2).^2+(1-x2).^2,x3=-2.048,y3=2.048,z3=100*(x3.*x3-y3).^2+(1-x3).^2,x4=-2.048,y4=-2.048,z4=100*(x4.*x4-y4).^2+(1-x4).^2%axis([-2.048,2.048,-2.048,2.048,-4000,4000]);。
遗传算法matlab程序编码>> function ret=Code(lenchrom,bound)%本函数将变量编码成染色体,用于随机初始化一个种群% lenchrom input : 染色体长度% bound input : 变量的取值范围% ret output: 染色体的编码值flag=0;while flag==0pick=rand(1,length(lenchrom));ret=bound(:,1)'+(bound(:,2)-bound(:,1))'.*pick; %线性插值,编码结果以实数向量存入ret 中flag=test(lenchrom,bound,ret); %检验染色体的可行性end交叉操作function ret=Cross(pcross,lenchrom,chrom,sizepop,bound)%本函数完成交叉操作% pcorss input : 交叉概率% lenchrom input : 染色体的长度% chrom input : 染色体群% sizepop input : 种群规模% ret output : 交叉后的染色体for i=1:sizepop %每一轮for循环中,可能会进行一次交叉操作,染色体是随机选择的,交叉位置也是随机选择的,%但该轮for循环中是否进行交叉操作则由交叉概率决定(continue 控制)% 随机选择两个染色体进行交叉pick=rand(1,2);while prod(pick)==0pick=rand(1,2);endindex=ceil(pick.*sizepop);% 交叉概率决定是否进行交叉pick=rand;while pick==0pick=rand;endif pick>pcrosscontinue;endflag=0;while flag==0% 随机选择交叉位pick=rand;while pick==0pick=rand;endpos=ceil(pick.*sum(lenchrom)); %随机选择进行交叉的位置,即选择第几个变量进行交叉,注意:两个染色体交叉的位置相同pick=rand; %交叉开始v1=chrom(index(1),pos);v2=chrom(index(2),pos);chrom(index(1),pos)=pick*v2+(1-pick)*v1;chrom(index(2),pos)=pick*v1+(1-pick)*v2; %交叉结束flag1=test(lenchrom,bound,chrom(index(1),:)); %检验染色体1的可行性flag2=test(lenchrom,bound,chrom(index(2),:)); %检验染色体2的可行性if flag1*flag2==0flag=0;else flag=1;end %如果两个染色体不是都可行,则重新交叉endendret=chrom;解码function ret=Decode(lenchrom,bound,code,opts)% 本函数对染色体进行解码% lenchrom input : 染色体长度% bound input : 变量取值范围% code input :编码值% opts input : 解码方法标签% ret output: 染色体的解码值switch optscase 'binary' % binary codingfor i=length(lenchrom):-1:1data(i)=bitand(code,2^lenchrom(i)-1); %并低十位,然后将低十位转换成十进制数存在data(i)里面code=(code-data(i))/(2^lenchrom(i)); %低十位清零,然后右移十位endret=bound(:,1)'+data./(2.^lenchrom-1).*(bound(:,2)-bound(:,1))'; %分段解码,以实数向量的形式存入ret中case 'grey' % grey codingfor i=sum(lenchrom):-1:2code=bitset(code,i-1,bitxor(bitget(code,i),bitget(code,i-1)));endfor i=length(lenchrom):-1:1data(i)=bitand(code,2^lenchrom(i)-1);code=(code-data(i))/(2^lenchrom(i));endret=bound(:,1)'+data./(2.^lenchrom-1).*(bound(:,2)-bound(:,1))'; %分段解码,以实数向量的形式存入ret中case 'float' % float codingret=code; %解码结果就是编码结果(实数向量),存入ret中end适应度函数function error = fun(x,inputnum,hiddennum,outputnum,net,inputn,outputn)%该函数用来计算适应度值%x input 个体%inputnum input 输入层节点数%outputnum input 隐含层节点数%net input 网络%inputn input 训练输入数据%outputn input 训练输出数据%error output 个体适应度值%提取w1=x(1:inputnum*hiddennum);B1=x(inputnum*hiddennum+1:inputnum*hiddennum+hiddennum);w2=x(inputnum*hiddennum+hiddennum+1:inputnum*hiddennum+hiddennum+hiddennum*out putnum);B2=x(inputnum*hiddennum+hiddennum+hiddennum*outputnum+1:inputnum*hiddennum+hid dennum+hiddennum*outputnum+outputnum);%网络进化参数net.trainParam.epochs=20;net.trainParam.lr=0.1;net.trainParam.goal=0.00001;net.trainParam.show=100;net.trainParam.showWindow=0;%网络权值赋值net.iw{1,1}=reshape(w1,hiddennum,inputnum);net.lw{2,1}=reshape(w2,outputnum,hiddennum);net.b{1}=reshape(B1,hiddennum,1);net.b{2}=B2;。
装配生产线任务平衡问题的遗传算法MATLAB源代码下面的源码实现了装配生产线任务平衡优化问题(ALB问题)的遗传算法,算法主要参考下面这篇文献,并对其进行了改进。
陈永卿,潘刚,李平.基于混合遗传算法的装配线平衡[J].机电工程,2008,25(4):60-62.。
function[BestX,BestY,BestZ,AllFarm,LC1,LC2,LC3,LC4,LC5]=GSAALB(M,N,Pm,Pd,K,t0,alpha,TaskP, TaskT,TaskV,RT,RV)% GreenSim团队——专业级算法设计&代写程序% 欢迎访问GreenSim团队主页→/greensim%% 装配生产线任务平衡问题的遗传算法%% 输入参数列表% M------------遗传算法进化代数% N------------种群规模,取偶数% Pm-----------变异概率调节参数% Pd-----------变异程度调节参数,0<Pd<1,越大,变异的基因位越多% K------------同一温度下状态跳转次数% T0-----------初始温度% Alpha--------降温系数% Beta---------浓度均衡系数% TaskP--------任务优先矩阵,n×n矩阵,Pij=1表示任务i需在j之前完成,Pij=0时任务i 和j没有优先关系% TaskT--------任务时间属性,n×1向量% TaskV--------任务体积属性,n×1向量% RT-----------时间节拍约束% RV-----------工位体积约束%% 输出参数列表% BestX--------最好个体的编码% BestY--------最好个体对应的装配方案% BestZ--------最好个体的目标函数值% LC1----------最优个体适应值的收敛曲线,M×1% LC2----------种群平均适应值的收敛曲线,M×1% LC3----------工位个数收敛曲线,M×1% LC4----------时间利用率及平衡度综合度量参数收敛曲线,M×1% LC5----------空间利用率及平衡度综合度量参数收敛曲线,M×1% AllFarm------各代种群的集合,M×1的细胞结构%% -----------------------初始化----------------------------------n=size(TaskP,1);[AA,BB]=QJHJ(TaskP);%调用子函数,建立每一个任务的前任务集和后任务集farm=Initialization(N,TaskP,AA,BB);%调用子函数,种群初始化%输出参数初始化BestX=zeros(1,n);BestY=zeros(1,n);BestZ=0;LC1=zeros(M,1);LC2=zeros(M,1);LC3=zeros(M,1);LC4=zeros(M,1);LC5=zeros(M,1);AllFarm=cell(M,1);%控制参数初始化m=1;%迭代计数器t=t0;%温度指示器BestPos=1;%初始时任意指定被保护个体%% -----------------------迭代过程---------------------------------while m<=M%设置停止条件%% ----------------------变异退火算子------------------------------for i=1:Nif rand>Pm&&i~=BestPos%如果随机数大于变异概率门限值,并且不属于保护个体,就对其实施变异I=farm(i,:);%取出该个体k=1;while k<=K%每一个温度下的状态转移次数%调用变异子函数J=Mutation(I,Pd,AA,BB);%调用计算适应值子函数[YI,ZI,FI,TGWI,VGWI,f1I,f2I]=Fitness(I,TaskT,TaskV,RT,RV);[YJ,ZJ,FJ,TGWJ,VGWJ,f1J,f2J]=Fitness(J,TaskT,TaskV,RT,RV);if FJ>FIfarm(i,:)=J;elseif rand<exp((FJ-FI)/(FI*t))farm(i,:)=J;elsefarm(i,:)=I;endk=k+1;endendend%% -----------------------交叉算子---------------------------------newfarm=zeros(size(farm));Ser=randperm(N);%用这个函数保证随机配对for i=1:2:(N-1)FA=farm(Ser(i),:);FB=farm(Ser(i+1),:);[SA,SB]=CrossOver(FA,FB);newfarm(i,:)=SA;newfarm(i+1,:)=SB;end%新旧种群合并FARM=[farm;newfarm];%% -----------------------选择复制--------------------------------- FIT_Y=zeros(2*N,n);FIT_Z=zeros(2*N,1);FIT_F=zeros(2*N,1);FIT_f1=zeros(2*N,1);FIT_f2=zeros(2*N,1);fit_Y=zeros(N,n);fit_Z=zeros(N,1);fit_F=zeros(N,1);fit_f1=zeros(N,1);fit_f2=zeros(N,1);for i=1:(2*N)XX=FARM(i,:);[Y,Z,F,TGW,VGW,f1,f2]=Fitness(XX,TaskT,TaskV,RT,RV);FIT_Y(i,:)=Y;FIT_Z(i)=Z;FIT_F(i)=F;FIT_f1(i)=f1;FIT_f2(i)=f2;endSer=randperm(2*N);for i=1:Nff1=FIT_F(Ser(2*i-1));ff2=FIT_F(Ser(2*i));if ff1>=ff2farm(i,:)=FARM(Ser(2*i-1),:);fit_Y(i,:)=FIT_Y(Ser(2*i-1),:);fit_Z(i)=FIT_Z(Ser(2*i-1));fit_F(i)=FIT_F(Ser(2*i-1));fit_f1(i)=FIT_f1(Ser(2*i-1));fit_f2(i)=FIT_f2(Ser(2*i-1));elsefarm(i,:)=FARM(Ser(2*i),:);fit_Y(i,:)=FIT_Y(Ser(2*i),:);fit_Z(i)=FIT_Z(Ser(2*i));fit_F(i)=FIT_F(Ser(2*i));fit_f1(i)=FIT_f1(Ser(2*i));fit_f2(i)=FIT_f2(Ser(2*i));endend%% -----------------------记录与更新------------------------------- maxF=max(fit_F);meanF=mean(fit_F);LC1(m)=maxF;LC2(m)=meanF;pos=find(fit_F==maxF);BestPos=pos(1);BestX=farm(BestPos,:);BestY=fit_Y(BestPos,:);BestZ=fit_Z(BestPos);LC3(m)=fit_Z(BestPos);LC4(m)=fit_f1(BestPos);LC5(m)=fit_f2(BestPos);AllFarm{m}=farm;disp(m);m=m+1;t=t*alpha;end源代码运行结果展示。
并行设计子任务调度的遗传算法MATLAB源代码在研究并行设计环境下设计子任务的调度问题时,针对该问题的描述所给出的假设往往各不相同不失一般性,我们假设一个大的总任务已经按照某种分解规则被分解为若干个子任务,这些子任务的求解存在有序性和可并行性,并且各子任务将由不同的设计单元来承担,而不同的设计单元对各子任务的求解效率是不尽相同的。
同时,假定需调度的子任务数多于可提供的设计单元数,并且各设计单元中的子任务之间的关系是非抢先式的,即在同一设计单元中的一个子任务没有完成之前,其他分配在该节点上的子任务不能开始执行。
优化目标是,在满足任务偏序图中各子任务间时序约束关系的条件下寻找一个任务调度策略,将子任务集中的n个子任务分别调度到设计单元集中的m个设计单元中去执行,并使得整个任务(总任务)的总的执行效率最高,即总的完成时间最短。
%%clcclearclose all%% 产生数据% GreenSim团队——专业级算法设计&代写程序% 欢迎访问GreenSim团队主页→/greensim%任务优先矩阵,n×n矩阵,Pij=1表示任务i需在j之前完成,Pij=0时任务i和j没有优先关系TaskP=zeros(8,8);TaskP(1,2)=1;TaskP(1,3)=1;TaskP(2,4)=1;TaskP(2,5)=1;TaskP(4,6)=1;TaskP(5,6)=1;TaskP(3,6)=1;TaskP(3,7)=1;TaskP(6,7)=1;TaskP(7,8)=1;%平均时间向量MeanTime=[20,30,85,15,20,105,52,40];%任务的效率矩阵EFF=[0.9, 0.8, 0.7, 0.9, 1.2, 0.9, 0.8, 1.2;1.0, 1.0, 1.1, 0.8, 0.9, 1.1, 1.0, 1.0;0.8, 1.1, 0.9, 1.1, 0.7, 1.2, 1.0, 0.9];[a,b]=size(EFF);Time=zeros(a,b);for i=1:aTime(i,:)=MeanTime./EFF(i,:);end[AA,BB]=QJHJ(TaskP);%%Iteration=1000;MinT=Inf;AllMinT=zeros(1,Iteration);Best_f1=zeros(1,b);Best_f2=zeros(1,b);Best_ST=zeros(1,b);Best_ET=zeros(1,b);for i=1:Iteration[f1,f2]=Initialization(1,TaskP,EFF,AA,BB);[ST,ET]=DiaoDu(f1,f2,AA,Time);T=max(ET);if T<MinTMinT=T;Best_f1=f1;Best_f2=f2;Best_ST=ST;Best_ET=ET;endAllMinT(i)=MinT;end%%figureDRAW(Best_f1,Best_f2,Best_ST,Best_ET);figureplot(AllMinT);function [ST,ET]=DiaoDu(f1,f2,AA,Time)%% 输入参数列表% GreenSim团队——专业级算法设计&代写程序% 欢迎访问GreenSim团队主页→/greensim % f1 用于描述任务的排序(符合约束)% f2 用于描述执行该任务的小组(机器)编号% AA 各任务的前任务集% Time 时间矩阵,各任务在各机器上的完成时间%% 输出参数列表% ST 各任务的开始时刻% ET 各任务的结束时刻%%[a,b]=size(Time);%a为分组个数,b为任务个数ET=zeros(1,b);ST=zeros(1,b);task=f1(1);%任务编号group=f2(1);%机器编号time=Time(group,task);ST(task)=0;ET(task)=time;for i=2:btask=f1(i);group=f2(i);time=Time(group,task);A T=AA{task};L=length(AT);et1=0;if L>0for j=1:Lif ET(A T(j))>et1et1=ET(AT(j));endendendP1=find(f2==group);BT=f1(P1);L2=length(BT);et2=0;if L2>0for j=1:L2if ET(BT(j))>et2et2=ET(BT(j));endendendST(task)=max([et1,et2]);ET(task)=ST(task)+time;end源代码运行结果展示。
流水线车间生产调度的遗传算法MATLAB源代码 n个任务在流水线上进行m个阶段的加工,每一阶段至少有一台机器且至少有一个阶段存在多台机器,并且同一阶段上各机器的处理性能相同,在每一阶段各任务均要完成一道工序,各任务的每道工序可以在相应阶段上的任意一台机器上加工,已知任务各道工序的处理时间,要求确定所有任务的排序以及每一阶段上机器的分配情况,使得调度指标(一般求Makespan)最小。
function [Zp,Y1p,Y2p,Y3p,Xp,LC1,LC2]=JSPGA(M,N,Pm,T,P) %-------------------------------------------------------------------------- % JSPGA.m % 流水线型车间作业调度遗传算法 % GreenSim团队——专业级算法设计&代写程序 % 欢迎访问GreenSim团队主页→blog.sina..cn/greensim %-------------------------------------------------------------------------- % 输入参数列表 % M 遗传进化迭代次数 % N 种群规模(取偶数) % Pm 变异概率 % T m×n的矩阵,存储m个工件n个工序的加工时间 % P 1×n的向量,n个工序中,每一个工序所具有的机床数目 % 输出参数列表 % Zp 最优的Makespan值 % Y1p 最优方案中,各工件各工序的开始时刻,可根据它绘出甘特图 % Y2p 最优方案中,各工件各工序的结束时刻,可根据它绘出甘特图 % Y3p 最优方案中,各工件各工序使用的机器编号 % Xp 最优决策变量的值,决策变量是一个实数编码的m×n矩阵 % LC1 收敛曲线1,各代最优个体适应值的记录 % LC2 收敛曲线2,各代群体平均适应值的记录 % 最后,程序还将绘出三副图片:两条收敛曲线图和甘特图(各工件的调度时序图) %第一步:变量初始化 [m,n]=size(T);%m是总工件数,n是总工序数 Xp=zeros(m,n);%最优决策变量 LC1=zeros(1,M);%收敛曲线1 LC2=zeros(1,N);%收敛曲线2 %第二步:随机产生初始种群 farm=cell(1,N);%采用细胞结构存储种群 for k=1:N X=zeros(m,n); for j=1:n for i=1:m X(i,j)=1+(P(j)-eps)*rand; end end farm{k}=X; end counter=0;%设置迭代计数器 while counter%第三步:交叉 newfarm=cell(1,N);%交叉产生的新种群存在其中 Ser=randperm(N); for i=1:2:(N-1) A=farm{Ser(i)};%父代个体 Manner=unidrnd(2);%随机选择交叉方式 if Manner==1 cp=unidrnd(m-1);%随机选择交叉点 %双亲双子单点交叉 a=[A(1:cp,:);B((cp+1):m,:)];%子代个体 b=[B(1:cp,:);A((cp+1):m,:)]; else cp=unidrnd(n-1);%随机选择交叉点 b=[B(:,1:cp),A(:,(cp+1):n)]; end newfarm{i}=a;%交叉后的子代存入newfarm newfarm{i+1}=b; end %新旧种群合并 FARM=[farm,newfarm]; %第四步:选择复制 FITNESS=zeros(1,2*N); fitness=zeros(1,N); plotif=0; for i=1:(2*N) X=FARM{i}; Z=COST(X,T,P,plotif);%调用计算费用的子函数 FITNESS(i)=Z; end %选择复制采取两两随机配对竞争的方式,具有保留最优个体的能力 Ser=randperm(2*N); for i=1:N f2=FITNESS(Ser(2*i)); if f1<=f2 farm{i}=FARM{Ser(2*i-1)}; fitness(i)=FITNESS(Ser(2*i-1)); else farm{i}=FARM{Ser(2*i)}; end end %记录最佳个体和收敛曲线 minfitness=min(fitness) meanfitness=mean(fitness) LC1(counter+1)=minfitness;%收敛曲线1,各代最优个体适应值的记录 LC2(counter+1)=meanfitness;%收敛曲线2,各代群体平均适应值的记录 pos=find(fitness==minfitness); Xp=farm{pos(1)}; %第五步:变异 for i=1:N if Pm>rand;%变异概率为Pm X=farm{i}; I=unidrnd(m); J=unidrnd(n); X(I,J)=1+(P(J)-eps)*rand; farm{i}=X; end end farm{pos(1)}=Xp; counter=counter+1 end %输出结果并绘图 figure(1); plotif=1; X=Xp; [Zp,Y1p,Y2p,Y3p]=COST(X,T,P,plotif); figure(2); plot(LC1); figure(3); plot(LC2); function [Zp,Y1p,Y2p,Y3p]=COST(X,T,P,plotif) % JSPGA的联子函数,用于求调度方案的Makespan值 % 输入参数列表 % X 调度方案的编码矩阵,是一个实数编码的m×n矩阵 % T m×n的矩阵,存储m个工件n个工序的加工时间 % P 1×n的向量,n个工序中,每一个工序所具有的机床数目 % plotif 是否绘甘特图的控制参数 % 输出参数列表 % Zp 最优的Makespan值 % Y1p 最优方案中,各工件各工序的开始时刻 % Y2p 最优方案中,各工件各工序的结束时刻 % Y3p 最优方案中,各工件各工序使用的机器编号 %第一步:变量初始化 [m,n]=size(X); Y1p=zeros(m,n); Y2p=zeros(m,n); Y3p=zeros(m,n); %第二步:计算第一道工序的安排 Q1=zeros(m,1); Q2=zeros(m,1); R=X(:,1);%取出第一道工序 Q3=floor(R);%向下取整即得到各工件在第一道工序使用的机器的编号 %下面计算各工件第一道工序的开始时刻和结束时刻 for i=1:P(1)%取出机器编号 pos=find(Q3==i);%取出使用编号为i的机器为其加工的工件的编号 lenpos=length(pos); if lenpos>=1 Q1(pos(1))=0; if lenpos>=2 for j=2:lenpos Q1(pos(j))=Q2(pos(j-1)); Q2(pos(j))=Q2(pos(j-1))+T(pos(j),1); end end end end Y1p(:,1)=Q1; Y3p(:,1)=Q3; %第三步:计算剩余工序的安排 for k=2:n R=X(:,k);%取出第k道工序 Q3=floor(R);%向下取整即得到各工件在第k道工序使用的机器的编号 %下面计算各工件第k道工序的开始时刻和结束时刻 for i=1:P(k)%取出机器编号 pos=find(Q3==i);%取出使用编号为i的机器为其加工的工件的编号 lenpos=length(pos); if lenpos>=1 EndTime=Y2p(pos,k-1);%取出这些机器在上一个工序中的结束时刻 POS=zeros(1,lenpos);%上一个工序完成时间由早到晚的排序 for jj=1:lenpos POS(jj)=ppp(1); EndTime(ppp(1))=Inf; end %根据上一个工序完成时刻的早晚,计算各工件第k道工序的开始时刻和结束时刻 Q1(pos(POS(1)))=Y2p(pos(POS(1)),k-1); Q2(pos(POS(1)))=Q1(pos(POS(1)))+T(pos(POS(1)),k);%前一个工件的结束时刻 if lenpos>=2 for j=2:lenpos Q1(pos(POS(j)))=Y2p(pos(POS(j)),k-1);%预定的开始时刻为上一个工序的结束时刻 if Q1(pos(POS(j)))还早 Q1(pos(POS(j)))=Q2(pos(POS(j-1))); end end end end end Y1p(:,k)=Q1; Y2p(:,k)=Q2; Y3p(:,k)=Q3; end %第四步:计算最优的Makespan值 Y2m=Y2p(:,n); Zp=max(Y2m); %第五步:绘甘特图 if plotif for i=1:m for j=1:n mPoint1=Y1p(i,j); mPoint2=Y2p(i,j); mText=m+1-i; PlotRec(mPoint1,mPoint2,mText); Word=num2str(Y3p(i,j)); %text(0.5*mPoint1+0.5*mPoint2,mText-0.5,Word); hold on x1=mPoint1;y1=mText-1; x2=mPoint2;y2=mText-1; x4=mPoint1;y4=mText; %fill([x1,x2,x3,x4],[y1,y2,y3,y4],'r'); fill([x1,x2,x3,x4],[y1,y2,y3,y4],[1,0.5,1]); text(0.5*mPoint1+0.5*mPoint2,mText-0.5,Word); end end end function PlotRec(mPoint1,mPoint2,mText)