1.rar_Benders分解法_Volterra_fredholm_三重积分_可靠度优化

function [chrx,MOC,PkG,LayoutG] = GA4WorkshopDeviceLayout %chrx:存储每一代的最优解;MOC:存储每一代的最优解的目标函数值, %PkG存储每一代的最优解的适应度， %LayoutG,存储每一代最优设备排列方式 DeviceSize = [1.2 1.0;0.5 0.4;1.7 0.8;0.7 0.7;1.0 0.8;0.8 0.8;... 0.5 0.8;1.4 0.8;0.8 0.8;1.2 0.6]; Pij = [0 8 2 1 6 4 3 2 3 5;8 0 3 2 6 7 3 4 2 5;2 3 0 3 4 2 4 3 1 6; 1 2 3 0 4 3 2 7 2 2;6 6 4 4 0 5 2 3 4 1;4 7 2 3 5 0 5 6 3 2; 3 3 4 2 2 5 0 4 4 4;2 4 3 7 3 6 4 0 3 2;3 2 1 2 4 3 4 3 0 2;... 5 5 6 2 1 2 4 2 2 0]; Qij = [0 3 2 1 2 1 1 0 2 1;3 0 2 0 1 1 3 2 0 1;2 2 0 1 1 0 1 3 2 1; 1 0 1 0 2 2 0 3 0 1;2 1 1 2 0 1 1 1 1 0;1 1 0 2 1 0 2 1 3 2; 1 3 1 0 1 2 0 2 3 3;0 2 3 3 1 1 2 0 1 3;2 0 2 0 1 3 3 1 0 2;... 1 1 1 1 0 2 3 3 2 0]; hi0 = [1 2 1.5 1.5 1.3 1.2 1 1.8 2 1]; hij = [0 1 1.8 1.8 1.2 1.4 2 1.3 1.2 1;1 0 2 2 1 1.5 1.7 1.6 1.5 1; 1.8 2 0 1 1 2 1 1.5 1.3 1;1.8 2 1 0 1 2 1.8 1.8 1 1; 1.2 1 1 1 0 1.7 1 1.7 2 1;1.4 1.5 2 2 1.7 0 2 1 1 1; 2 1.7 1 1.8 1 2 0 1.2 1 1;1.3 1.6 1.5 1.8 1.7 1 1.2 0 1.8 1; 1.2 1.5 1.3 1 2 1 1 1.8 0 1;1 1 1 1 1 1 1 1 1 0]; s = 2; s0 = 1.5; nchr = 50;%种群数量 G = 200;%代终止数 pc = 0.6;%交叉概率 pm = 0.1;%变异概率 T = 500; H = 8;%车间宽 r = 10; Umin = 0; Umax = 1.5; %% 计算距离矩阵的函数(采用nested function结构，以function+end标识,位置可以出 %现在主函数体任何位置) % chr：染色体。Dij：返回的距离矩阵 function [Dij,nRow,Layout] = GenerateDistanceMatrix(chr) m = chr(1:10); delta = chr(11:20); %记录每行排列的设备号，由于行数不定，但最多不可能超过10行，所以用一个10*1的 %cell数组记录。 Layout = cell(10,1); k = 1;%染色体设备序列m中第一个设备号 nRow = 1;%排列的行数 while k<10 %while循环：按照自动换行策略判断设备是否排列完 for kk=k:10 if sum(DeviceSize(k:kk,1)) + sum(delta(k:kk)) + ... sum( hij( sub2ind( size(hij), m(k:kk-1), m(k+1:kk) )))... + sum( hi0( m([k,kk]) ) ) > 10 %判断是否该换行，第一项：设备总长，第二项：净间距和，第三项和第 %四项相应的hij和hi0之和 Layout{nRow} = m(k:kk-1); k = kk; break; end end nRow = nRow+1; %判断是否到最后一行，并布置最后一行的设备 if kk == 10 Layout{nRow} = m(kk); elseif sum(DeviceSize(kk:end,1)) + sum(delta(kk:end)) +... sum( hij( sub2ind( size(hij), m(kk:end-1), m(kk+1:end) )))... + sum( hi0( m([kk,end]) ) ) <= 10 Layout{nRow} = m(kk:end); break end end x = zeros(10,1);%存储各设备x坐标 y = zeros(10,1);%存储各设备y坐标 for nr = 1:nRow x( Layout{nr}(1) ) = hi0( Layout{nr}(1) ) + delta( Layout{nr}(1) ); y(Layout{nr}) = (nr-1)*s+s0; if length(Layout{nr}) == 1 %某行只有一个设备，上面已经得到设备的横纵坐标，直接循环下一个nr continue; end for nc = 2:length(Layout{nr}) x(Layout{nr}(nc)) = x(Layout{nr}(nc-1)) +... mean( DeviceSize(Layout{nr}([nc-1,nc]),1) ) + ... hij( Layout{nr}(nc-1),Layout{nr}(nc) )+delta(Layout{nr}(nc)); end end Dij = zeros(10); for ii=1:10 for jj=1:10 Dij(ii,jj) = abs(x(ii) - x(jj))+abs( y(ii)-y(jj) ); end end end %% 随机生成nchr个初始种群,用已有设备序列替换掉前三个染色体的设备序列段 %M,由于是对随机数矩阵每行排序后新矩阵中的数字在原矩阵的序号，所以M相当于随机生成了nchr %个1到10的随机排序 [sorted,M] = sort(rand(nchr,10),2); CHR = [M,unifrnd(Umin,Umax,nchr,10)] ; for ki = 1:3 %用给的设备序列替换掉前三个染色体的设备序列段 CHR(ki,1:10) = [2 8 5 4 1 7 3 9 10 6]; end %% 目标函数以及适应度函数，设备排列，根据适应度函数特点，将目标函数和适应度函数 %以及设备排列用一个函数实现，减少重复计算 function [C,Pk,Layout] = ObjfunEvalfun(chr) [Dij,nRow,Layout] = GenerateDistanceMatrix(chr); C = sum( sum(Pij.*Qij.*Dij) );%目标函数值 Pk = 1/( C+T*( (nRow-1)*s+s0>H ) );%适应度函数值 end %% 进化G次,PkG存储每一代的最优解的适应度，MOC,存储每一代的最优解的目标函数值， %chrx存储每一代的最优解，C，进化过程中每代染色体的目标函数值，Pk，进化过程中每代染色体 %的适应度，LayoutG,存储每一代最优设备排列方式，Layoutnchr,每一代中种群各个染色体的 %设备排列方式 PkG = zeros(G,1); MOC = zeros(G,1); chrx=zeros(G,20); LayoutG = cell(G,1); Layoutnchr = cell(nchr,1); C = zeros(nchr,1); Pk = zeros(nchr,1); for g = 1:G %% 计算目标函数值以及适应度 for kn = 1:nchr [C(kn),Pk(kn),Layoutnchr{kn}] = ObjfunEvalfun(CHR(kn,:)); end [sortPk, ind] = sort(Pk); %进化的代数大于1，但是最大适应度值不如上一代的优。用上一代的最优结果替代本代最优结果. if g>1 && sortPk(end) < PkG(g-1) PkG(g) = PkG(g-1); chrx(g,:) = chrx(g-1,:); MOC(g) = MOC(g-1); LayoutG{g} = LayoutG{g-1}; else PkG(g) = sortPk(end); chrx(g,:) = CHR(ind(end),:); MOC(g) = C(ind(end)); LayoutG{g} = Layoutnchr{ind(end)}; end %% 选择操作 SelectP = Pk/sum(Pk);%选择概率 %01区间划分（各区间长度依次等于各个体选择的概率比例） Interval01Divide=cumsum(SelectP); Test=rand(nchr,1);%生成nchr（50）个0，1随机数 %随机数落入的区间序号，表示该区间代表的个体被选中 TestResult=arrayfun(@(x) find(Interval01Divide>=x,1,'first'),Test); Frequency=accumarray(TestResult,1);%各个体选中的次数，进而决定各个体复制的次数 Chr=cell(length(Frequency),1); for tt=1:length(Frequency) %根据Frequency每个分量决定相应个体的复制次数 Chr{tt}=repmat(CHR(tt,:),Frequency(tt),1); end Chr=cell2mat(Chr); %% 该部分是交叉操作 RandMatch=randperm(nchr);%随机配对 for tt=1:nchr/2 if rand>1-pc %满足交叉概率,下面交叉 a1 = Chr(RandMatch(2*tt-1),1:10);%a1,a2进行交叉的染色体设备序列段 a2 = Chr(RandMatch(2*tt),1:10); b1 = Chr(RandMatch(2*tt-1),11:20);%b1,b2进行交叉的染色体delta段 b2 = Chr(RandMatch(2*tt),11:20); Points = sort(unidrnd(10,1,2)); %交叉a1,a2 for uu = 0:range(Points) a1(a1==a2(Points(1)+uu)) = a1(Points(1)+uu); temp = a1(Points(1)+uu); a1(Points(1)+uu) = a2(Points(1)+uu); a2(a2==temp) = a2(Points(1)+uu); a2(Points(1)+uu) = temp; end %交叉b1,b2 alpha = unidrnd(2,1,10)-1; temp2 = b1.*alpha+b2.*(1-alpha); b2 = b1.*(1-alpha)+b2.*alpha; b1 = temp2; Chr(RandMatch(2*tt-1),:) = [a1,b1]; Chr(RandMatch(2*tt),:) = [a2,b2]; end end %% 该部分是变异操作 for tt = 1:nchr if rand<pm %满足变异概率 pos = unidrnd(10,1)+10; %生成delta变异位置(10个位置中的随机一个) delta_r = unifrnd(Umin,Umax,1,r); %随机生成r个delta PkTemp = zeros(1,r); for tr = 1:r Chr(tt,pos) = delta_r(tr); [Ctemp,PkTemp(tr)] = ObjfunEvalfun(Chr(tt,:)); end [maxPktemp,indmax] = max(PkTemp); Chr(tt,pos) = delta_r(indmax); end end %% 种群更新 CHR = Chr; end plot(1:G,MOC); xlabel('进化次数'); ylabel('最优目标函数值'); title('进化次数同最优函数值曲线') end