天天看点

2021数学建模国赛C题程序 跑通版 程序开源

原文链接:

2021数学建模国赛C题思路 生产企业原材料的订购与运输 第一版思路 思路开源 已修订 程序_您好啊数模君的博客-CSDN博客这里有百种算法出处整理,本题算法可从上面找取:给裸赛的家人们整理了百种算法出处https://mp.weixin.qq.com/s/OhWRCeep885MuyhMhvdiOw这道题是优化问题,先看看题目,以下为题目告知的条件和目标函数条件:48周;24周的原材料订购和转运计划,就需要达到48周的用料企业每周的产能为2.82万立方米,每立方米产品需消耗A类原材料0.6立方米,或B类原材料0.66立方米,或C类原材料0.72立方米;企业对供应商实际提供的原材料总是全部收购;

2021数学建模国赛C题程序 跑通版 程序开源

https://blog.csdn.net/qq_39899679/article/details/120212554

程序文件版下载链接

链接:https://pan.baidu.com/s/13aSy2-hlkLa7Ps8wknYY1Q 

提取码:gap4 

第一问

clear
X1=xlsread('附件1 近5年402家供应商的相关数据.xlsx','企业的订货量(m³)');
X2=xlsread('附件1 近5年402家供应商的相关数据.xlsx','供应商的供货量(m³)');
[~,X3,~]=xlsread('附件1 近5年402家供应商的相关数据.xlsx','企业的订货量(m³)');
X3=string(X3);
X3=X3(2:end,2);
X3(find(X3=="A"))=1.3889;
X3(find(X3=="B"))=1.3774;
X3(find(X3=="C"))=1.3889;
%从供应角度看供应商的供应情况

Y(:,1)=1./std(X2')';%供应商供应的稳定性,对历史供应数据,做方差,负向指标,可以做个倒数
Y(:,2)=double(X3);%供应商材料经济效益,正向指标
Y(:,3)=sum(X2,2)./sum(X1,2);%历史供应量/历史订单量比,正向指标
%从企业角度看对供应商的认可情况
Y(:,4)=max(X1')';%供应商最大订单量,正向指标
Y(:,5)=mean(X1,2);%供应商平均订单量,正向指标
Y(:,6)=1./std(X1')';%供应商订单量的稳定性,对历史供应数据,做方差,负向指标,可以做个倒数
%还有其他指标可自加
Y=log(1+Y);%对指标ln函数映射,以减少数据本身差距,将评分重心转移至权重
Y(:,1)=log(100+Y(:,1));%对指标ln函数映射,以减少数据本身差距,将评分重心转移至权重
Y=mapminmax(Y',0.1,1)';
%% 熵权法
[n,m]=size(Y);
for i=1:n
    for j=1:m
        p(i,j)=Y(i,j)/sum(Y(:,j));
    end
end
%% 计算第 j 个指标的熵值 e(j)
k=1/log(n);
for j=1:m
    e(j)=-k*sum(p(:,j).*log(p(:,j)));
end
d=ones(1,m)-e; % 计算信息熵冗余度
w=d./sum(d); % 求权值 w
f=Y*w';
[f,paixu]=sort(f,'descend');
aa=paixu(1:50);
save aa aa
           

第二问

2021数学建模国赛C题程序 跑通版 程序开源
clear
%% 第一部分请自行修改为第一问的Top50供应商
%直接运行程序十有八九会报错,因为是随机选的50,请将aa矩阵更改为你们第一问的结果Top50的编号
[~,X1,~]=xlsread('附件1 近5年402家供应商的相关数据.xlsx','企业的订货量(m³)');
X1=string(X1);
X1=X1(2:end,2);
X2=xlsread('附件1 近5年402家供应商的相关数据.xlsx','供应商的供货量(m³)');
X3=xlsread('附件2 近5年8家转运商的相关数据.xlsx');
X3=mean(X3,2);%这里我直接取平均损耗了
XXX3=X3;
[X3,x3]=sort(X3);
X3=[X3,6000.*ones(8,1)];%加上最大转运量,方便调用
Y=zeros(size(X1,1),3);%记录单价、单位原料转化为产能的比例、供应商最大供应能力
Y(find(X1=="A"),1)=1.2;
Y(find(X1=="B"),1)=1.1;
Y(find(X1=="C"),1)=1;
Y(find(X1=="A"),2)=1/0.6;
Y(find(X1=="B"),2)=1/0.66;
Y(find(X1=="C"),2)=1/0.72;
Y(:,3)=max(X2')';%供应商最大供应能力,正向指标
%我这里随机选择50个,你们的话就用第一问结果
clear X2
load aa
X1=X1(aa,:);Y=Y(aa,:);

%% 寻优部分
num=100;%个体数
T0=100;%初始温度
q=0.95;%降温系数
Tmin=1;%最低温度
 
%初始化种群
x=[];
T=[];
F=[];
for i=1:num
 
    YL=48*28200; %总产能
    tt1=[];
    tt2=[];
    while YL>0 %直到满足产能为止
 
    if YL>56400
        z=56400;
    else
        z=YL;
    end
    t1=zeros(50,1);
    %这里是随机给个初始值,方便直接进入while循环
    k=[];k=randi([3,50]);%最多选择50个供应商
    a=[];a=randperm(50);a=a(1:k);%随机选取供应商
    a=sort(a);
    while sum(Y(a,2).*Y(a,3).*(1-max(X3(:,1))/100))<z %判断其最大供应量转化为产能是否达标,不能则重新随机
        k=[];k=randi([3,50]);%最多选择50个供应商
        a=[];a=randperm(50);a=a(1:k);%随机选取供应商
        a=sort(a);
    end
    
    Min=0;%记录供应量
    m=0;%记录对应产能
    n=0;
    b=[];[~,b]=sort(X1(a),'descend');b=a(b);%按CBA顺序排序,确定Min
    while m<z
        n=n+1;
        if m+Y(b(n),2)*Y(b(n),3)<z
            m=m+Y(b(n),2)*Y(b(n),3)*(1-max(X3(:,1))/100);
            Min=Min+Y(b(n),3);
        else
            Min=Min+(z-m)/Y(b(n),2);
            m=z;
        end
    end
    c=[];c=Min/sum(Y(a,3));

    d=50000;
    while sum(d)>48000
        d=round((c+(1-c).*rand(k,1)).*Y(a,3));%由于从不同供应商那里订购同一材料价格相同,实则最后就是ABC原料的订购比
    end
    t1(a,1)=d;%订单
    %转运,根据思路将A优先安排给损耗最小的转运商
    bb=[];[~,bb]=sort(X1(a));bb=a(bb);
    XX3=X3(:,2);
    t2=zeros(50,8);
    for j=1:length(bb)
        f=[];f=t1(bb(j));
        while f>0
            e=[];e=find(XX3>0);
            if XX3(e(1))>=f
                t2(bb(j),x3(e(1)))=f;
                XX3(e(1))=XX3(e(1))-f;
                f=0;
            else
                t2(bb(j),x3(e(1)))=XX3(e(1));
                f=f-XX3(e(1));
                XX3(e(1))=0;
            end
        end
    end
    YL=YL-sum(t2.*Y(:,2)*(1-XXX3./100)); %将一周的供应量换算为产能
    tt1=[tt1,t1];%格式按订购方案表格排列
    tt2=[tt2,t2];%格式按转运方案表格排列
    end
    %记录变量及目标函数
    x(i,1)=k;
    T{i,1}=tt1;
    T{i,2}=tt2;
    F(i,1)=k;%供应商数
    F(i,2)=2*sum(sum(tt2));%总费用
end
[TT,x,T]=ns2_1(x,T,F(:,1),F(:,2));
bestf=TT(1,:);
bestx=x(1,:);
bestT=T(1,:);
figure
plot(TT(:,1),TT(:,2),'b*')
 
while T0>Tmin
    %初始化种群
    xt=[];
    Tt=[];
    Ft=[];
    for i=1:num
        
        YL=48*28200;
        tt1=[];
        tt2=[];
        while YL>0
        if YL>56400
            z=56400;
        else
            z=YL;
        end
        t1=zeros(50,1);
        k=[];k=randi([3,50]);%这里是随机给个初始值,方便直接进入while循环
        a=[];a=randperm(50);a=a(1:k);
        a=sort(a);
        while sum(Y(a,2).*Y(a,3).*(1-max(X3(:,1))/100))<z %判断其最大供应量转化为产能是否达标,不能则重新随机
            k=[];k=randi([3,50]);
            a=[];a=randperm(50);a=a(1:k);
            a=sort(a);
        end
    
        Min=0;%记录供应量
        m=0;%记录对应产能
        n=0;
        b=[];[~,b]=sort(X1(a),'descend');b=a(b);
        while m<z
            n=n+1;
            if m+Y(b(n),2)*Y(b(n),3)<z
                m=m+Y(b(n),2)*Y(b(n),3)*(1-max(X3(:,1))/100);
                Min=Min+Y(b(n),3);
            else
                Min=Min+(z-m)/Y(b(n),2);
                m=z;
            end
        end
        c=[];c=Min/sum(Y(a,3));
        
        d=50000;
        while sum(d)>48000
            d=round((c+(1-c).*rand(k,1)).*Y(a,3));%由于从不同供应商那里订购同一材料价格相同,实则最后就是ABC原料的订购比
        end
        t1(a,1)=d;%订单
        %转运
        bb=[];[~,bb]=sort(X1(a));bb=a(bb);
        XX3=X3(:,2);
        t2=zeros(50,8);
        for j=1:length(bb)
            f=[];f=t1(bb(j));
            while f>0
                e=[];e=find(XX3>0);
                if XX3(e(1))>=f
                    t2(bb(j),x3(e(1)))=f;
                    XX3(e(1))=XX3(e(1))-f;
                    f=0;
                else
                    t2(bb(j),x3(e(1)))=XX3(e(1));
                    f=f-XX3(e(1));
                    XX3(e(1))=0;
                end
            end
        end
        YL=YL-sum(t2.*Y(:,2)*(1-XXX3./100));
        tt1=[tt1,t1];
        tt2=[tt2,t2];
        end
        %记录变量及目标函数
        xt(i,1)=k;
        Tt{i,1}=tt1;
        Tt{i,2}=tt2;
        Ft(i,1)=k;%供应商数
        Ft(i,2)=2*sum(sum(tt2));
    end
    [TTt,xt,Tt]=ns2_1([x;xt],[T;Tt],[F(:,1);Ft(:,1)],[F(:,2);Ft(:,2)]);
    TTt=TTt(1:num,:);
    xt=xt(1:num,:);
    Tt=Tt(1:num,:);
    for kk=1:num
        delta=TTt(kk,1)*TTt(kk,2)-TT(kk,1)*TT(kk,2);
        if delta<0
            TT(kk,:)=TTt(kk,:);
            x(kk,:)=xt(kk,:);
            T(kk,:)=Tt(kk,:);
        else
            P=exp(-delta/T0);
            if P>rand
                TT(kk,:)=TTt(kk,:);
                x(kk,:)=xt(kk,:);
                T(kk,:)=Tt(kk,:);
            end
        end
    end
    bestf=TT(1,:);
    bestx=x(1,:);
    bestT=T(1,:);
    T0=T0*q;
end
hold on
plot(TT(:,1),TT(:,2),'r*')
legend('初始分布','最优分布')
title('模拟退火寻优-解分布')
xlabel('F1 供应商数')
ylabel('F2 总费用')
%% 最后的结果bestf为最优方案的目标函数,bestx为最少供应商数,bestT中有两个矩阵,均以按附件AB表格排好格式,直接复制过去即可
           

第三问

2021数学建模国赛C题程序 跑通版 程序开源
clear
%% 第一部分请自行修改为第一问的Top50供应商
%直接运行程序十有八九会报错,因为是随机选的50,请将aa矩阵更改为你们第一问的结果Top50的编号
[~,X1,~]=xlsread('附件1 近5年402家供应商的相关数据.xlsx','企业的订货量(m³)');
X1=string(X1);
X1=X1(2:end,2);
X2=xlsread('附件1 近5年402家供应商的相关数据.xlsx','供应商的供货量(m³)');
X3=xlsread('附件2 近5年8家转运商的相关数据.xlsx');
X3=mean(X3,2);%这里我直接取平均损耗了
XXX3=X3;
[X3,x3]=sort(X3);
X3=[X3,6000.*ones(8,1)];%加上最大转运量,方便调用
Y=zeros(size(X1,1),3);%记录单价、单位原料转化为产能的比例、供应商最大供应能力
Y(find(X1=="A"),1)=1.2;
Y(find(X1=="B"),1)=1.1;
Y(find(X1=="C"),1)=1;
Y(find(X1=="A"),2)=1/0.6;
Y(find(X1=="B"),2)=1/0.66;
Y(find(X1=="C"),2)=1/0.72;
Y(:,3)=max(X2')';%供应商最大供应能力,正向指标
%我这里随机选择50个,你们的话就用第一问结果
clear X2
load aa
X1=X1(aa,:);Y=Y(aa,:);

%% 寻优部分
num=300;%个体数
T0=100;%初始温度
q=0.95;%降温系数
Tmin=1;%最低温度

%初始化种群
x=[];
T=[];
F=[];
for i=1:num

    YL=48*28200; %总产能
    tt1=[];
    tt2=[];
    while YL>0 %直到满足产能为止

    if YL>56400
        z=56400;
    else
        z=YL;
    end
    t1=zeros(50,1);
    %这里是随机给个初始值,方便直接进入while循环
    k=[];k=randi([3,50]);%最多选择50个供应商
    a=[];a=randperm(50);a=a(1:k);%随机选取供应商
    a=sort(a);
    while sum(Y(a,2).*Y(a,3).*(1-max(X3(:,1))/100))<z %判断其最大供应量转化为产能是否达标,不能则重新随机
        k=[];k=randi([3,50]);%最多选择50个供应商
        a=[];a=randperm(50);a=a(1:k);%随机选取供应商
        a=sort(a);
    end
    Min=0;%记录供应量
    m=0;%记录对应产能
    n=0;
    b=[];[~,b]=sort(X1(a),'descend');b=a(b);%按CBA顺序排序,确定Min
    while m<z
        n=n+1;
        if m+Y(b(n),2)*Y(b(n),3)<z
            m=m+Y(b(n),2)*Y(b(n),3)*(1-max(X3(:,1))/100);
            Min=Min+Y(b(n),3);
        else
            Min=Min+(z-m)/Y(b(n),2);
            m=z;
        end
    end
    c=[];c=Min/sum(Y(a,3));
    d=50000;
    while sum(d)>48000
        d=round((c+(1-c).*rand(k,1)).*Y(a,3));%由于从不同供应商那里订购同一材料价格相同,实则最后就是ABC原料的订购比
    end
    t1(a,1)=d;%订单
    %转运,根据思路将A优先安排给损耗最小的转运商
    bb=[];[~,bb]=sort(X1(a));bb=a(bb);
    XX3=X3(:,2);
    t2=zeros(50,8);
    for j=1:length(bb)
        f=[];f=t1(bb(j));
        while f>0
            e=[];e=find(XX3>0);
            if XX3(e(1))>=f
                t2(bb(j),x3(e(1)))=f;
                XX3(e(1))=XX3(e(1))-f;
                f=0;
            else
                t2(bb(j),x3(e(1)))=XX3(e(1));
                f=f-XX3(e(1));
                XX3(e(1))=0;
            end
        end
    end
    YL=YL-sum(t2.*Y(:,2)*(1-XXX3./100)); %将一周的供应量换算为产能
    tt1=[tt1,t1];%格式按订购方案表格排列
    tt2=[tt2,t2];%格式按转运方案表格排列
    end
    %记录变量及目标函数
    x(i,1)=k;
    T{i,1}=tt1;
    T{i,2}=tt2;
    F(i,1)=k;%供应商数
    F(i,2)=2*sum(sum(tt2));%总费用
    F(i,3)=sum(sum(tt1(find(X1=="A"),:)));%A材料
    F(i,4)=sum(sum(tt1(find(X1=="C"),:)));%C材料
end
[TT,x,T]=ns2_2(x,T,F(:,1),F(:,2),F(:,3),F(:,4));
bestf=TT(1,:);
bestx=x(1,:);
bestT=T(1,:);
figure
plot3(TT(:,1),TT(:,2),TT(:,3),'b*')

while T0>Tmin
    %初始化种群
    xt=[];
    Tt=[];
    Ft=[];
    for i=1:num
        
        YL=48*28200;
        tt1=[];
        tt2=[];
        while YL>0
        if YL>56400
            z=56400;
        else
            z=YL;
        end
        t1=zeros(50,1);
        k=[];k=randi([3,50]);
        a=[];a=randperm(50);a=a(1:k);
        a=sort(a);%这里是随机给个初始值,方便直接进入while循环
        while sum(Y(a,2).*Y(a,3).*(1-max(X3(:,1))/100))<z %判断其最大供应量转化为产能是否达标,不能则重新随机
            k=[];k=randi([3,50]);
            a=[];a=randperm(50);a=a(1:k);
            a=sort(a);
        end
        Min=0;%记录供应量
        m=0;%记录对应产能
        n=0;
        b=[];[~,b]=sort(X1(a),'descend');b=a(b);
        while m<z
            n=n+1;
            if m+Y(b(n),2)*Y(b(n),3)<z
                m=m+Y(b(n),2)*Y(b(n),3)*(1-max(X3(:,1))/100);
                Min=Min+Y(b(n),3);
            else
                Min=Min+(z-m)/Y(b(n),2);
                m=z;
            end
        end
        c=[];c=Min/sum(Y(a,3));
        d=50000;
        while sum(d)>48000
            d=round((c+(1-c).*rand(k,1)).*Y(a,3));%由于从不同供应商那里订购同一材料价格相同,实则最后就是ABC原料的订购比
        end
        t1(a,1)=d;%订单
        %转运
        bb=[];[~,bb]=sort(X1(a));bb=a(bb);
        XX3=X3(:,2);
        t2=zeros(50,8);
        for j=1:length(bb)
            f=[];f=t1(bb(j));
            while f>0
                e=[];e=find(XX3>0);
                if XX3(e(1))>=f
                    t2(bb(j),x3(e(1)))=f;
                    XX3(e(1))=XX3(e(1))-f;
                    f=0;
                else
                    t2(bb(j),x3(e(1)))=XX3(e(1));
                    f=f-XX3(e(1));
                    XX3(e(1))=0;
                end
            end
        end
        YL=YL-sum(t2.*Y(:,2)*(1-XXX3./100));
        tt1=[tt1,t1];
        tt2=[tt2,t2];
        end
        %记录变量及目标函数
        xt(i,1)=k;
        Tt{i,1}=tt1;
        Tt{i,2}=tt2;
        Ft(i,1)=k;%供应商数
        Ft(i,2)=2*sum(sum(tt2));%总费用
        Ft(i,3)=sum(sum(tt1(find(X1=="A"),:)));%A材料
        Ft(i,4)=sum(sum(tt1(find(X1=="C"),:)));%C材料
    end
    [TTt,xt,Tt]=ns2_2([x;xt],[T;Tt],[F(:,1);Ft(:,1)],[F(:,2);Ft(:,2)],[F(:,3);Ft(:,3)],[F(:,4);Ft(:,4)]);
    TTt=TTt(1:num,:);
    xt=xt(1:num,:);
    Tt=Tt(1:num,:);
    for kk=1:num
        delta=TTt(kk,3)/(TTt(kk,1)*TTt(kk,2)*TTt(kk,4))-TT(kk,3)/(TT(kk,1)*TT(kk,2)*TT(kk,4));
        if delta<0
            TT(kk,:)=TTt(kk,:);
            x(kk,:)=xt(kk,:);
            T(kk,:)=Tt(kk,:);
        else
            P=exp(-delta/T0);
            if P>rand
                TT(kk,:)=TTt(kk,:);
                x(kk,:)=xt(kk,:);
                T(kk,:)=Tt(kk,:);
            end
        end
    end
    bestf=TT(1,:);
    bestx=x(1,:);
    bestT=T(1,:);
    T0=T0*q;
end
hold on
plot3(TT(:,1),TT(:,2),TT(:,3),'r*')
legend('初始分布','最优分布')
title('模拟退火寻优-解分布')
xlabel('F1 供应商数')
ylabel('F2 总费用')
zlabel('F3 A材料供应数')
%% 最后的结果bestf为最优方案的目标函数,bestx为最少供应商数,bestT中有两个矩阵,均以按附件AB表格排好格式,直接复制过去即可
           

其他自定义函数,需创建调用

function [TT,chrom,cchrom]=ns2(x,T,F1,F2)
% 快速非支配排序
a = 0;
T1 = [];
T2 = [];
chrom=x;
cchrom=T;
chrom1 = [];
chrom2 = [];
cchrom1 = [];
cchrom2 = [];
while a == 0 %根据被支配数进行分级和排序
    M = [];
    for i = 1:length(F1)
        M(i,1) = length(find(F1<F1(i,1)))+length(find(F2<F2(i,1)));%目标函数最小化这里为<,最大化改成>
    end
    b1 = [];
    b2 = [];
    [b1,b2] = sort(M); %b1返回从小到大排序,b2返回原始序号
    if  length(chrom)>0 && b1(1) == 0   %无被支配数进入一级用T1矩阵保存
        T1 = [T1;F1(b2(1)),F2(b2(1))];
        chrom1 = [chrom1;chrom(b2(1),:)];
        cchrom1 = [cchrom1;cchrom(b2(1),:)];
        F1(b2(1)) = [];
        F2(b2(1)) = [];
        chrom(b2(1),:) = [];
        cchrom(b2(1),:) = [];
    else  %有被支配数进入二级用T2矩阵保存
        a = 1;
        T2 = [F1,F2];
        chrom2 = chrom;
        cchrom2 = cchrom;
    end
end
T2 = T2(b2,:);
chrom2 = chrom2(b2,:);
cchrom2 = cchrom2(b2,:);
if size(T1,1) > 2  %T1矩阵不用进行拥挤度调整排序,直接对T2进行排序调整即可
    y = yongji(T1);%拥挤度
    for i = 2:size(T1,1)
            if y(i-1) > y(i)
                T1(i-1:1:i,:) = T1(i:-1:i-1,:); %根据拥挤度调整排序,如果后者优于前者则反转顺序
                chrom1(i-1:1:i,:) = chrom1(i:-1:i-1,:);
                cchrom1(i-1:1:i,:) = cchrom1(i:-1:i-1,:);
            end
    end
end

if length(T2) > 0  %T1矩阵不用进行拥挤度调整排序,直接对T2进行排序调整即可
    y = yongji(T2);%拥挤度
    for i = 2:size(T2,1)
        if b1(i) == b1(i-1)
            if y(i-1) > y(i)
                T2(i-1:1:i,:) = T2(i:-1:i-1,:); %根据拥挤度调整排序,如果后者优于前者则反转顺序
                chrom2(i-1:1:i,:) = chrom2(i:-1:i-1,:);
                cchrom2(i-1:1:i,:) = cchrom2(i:-1:i-1,:);
            end
        end
    end
end
%排序重组
TT = [T1;T2];
chrom = [chrom1;chrom2];
cchrom = [cchrom1;cchrom2];
           
function [TT,chrom,cchrom]=ns2(x,T,F1,F2,F3,F4)
% 快速非支配排序
a = 0;
T1 = [];
T2 = [];
chrom=x;
cchrom=T;
chrom1 = [];
chrom2 = [];
cchrom1 = [];
cchrom2 = [];
while a == 0 %根据被支配数进行分级和排序
    M = [];
    for i = 1:length(F1)
        M(i,1) = length(find(F1<F1(i,1)))+length(find(F2<F2(i,1)))+length(find(F3>F3(i,1)))+length(find(F4<F4(i,1)));%目标函数最小化这里为<,最大化改成>
    end
    b1 = [];
    b2 = [];
    [b1,b2] = sort(M); %b1返回从小到大排序,b2返回原始序号
    if  length(chrom)>0 && b1(1) == 0   %无被支配数进入一级用T1矩阵保存
        T1 = [T1;F1(b2(1)),F2(b2(1)),F3(b2(1)),F4(b2(1))];
        chrom1 = [chrom1;chrom(b2(1),:)];
        cchrom1 = [cchrom1;cchrom(b2(1),:)];
        F1(b2(1)) = [];
        F2(b2(1)) = [];
        F3(b2(1)) = [];
        F4(b2(1)) = [];
        chrom(b2(1),:) = [];
        cchrom(b2(1),:) = [];
    else  %有被支配数进入二级用T2矩阵保存
        a = 1;
        T2 = [F1,F2,F3,F4];
        chrom2 = chrom;
        cchrom2 = cchrom;
    end
end
T2 = T2(b2,:);
chrom2 = chrom2(b2,:);
cchrom2 = cchrom2(b2,:);
if size(T1,1) > 2  %T1矩阵不用进行拥挤度调整排序,直接对T2进行排序调整即可
    y = yongji(T1);%拥挤度
    for i = 2:size(T1,1)
            if y(i-1) > y(i)
                T1(i-1:1:i,:) = T1(i:-1:i-1,:); %根据拥挤度调整排序,如果后者优于前者则反转顺序
                chrom1(i-1:1:i,:) = chrom1(i:-1:i-1,:);
                cchrom1(i-1:1:i,:) = cchrom1(i:-1:i-1,:);
            end
    end
end

if length(T2) > 0  %T1矩阵不用进行拥挤度调整排序,直接对T2进行排序调整即可
    y = yongji(T2);%拥挤度
    for i = 2:size(T2,1)
        if b1(i) == b1(i-1)
            if y(i-1) > y(i)
                T2(i-1:1:i,:) = T2(i:-1:i-1,:); %根据拥挤度调整排序,如果后者优于前者则反转顺序
                chrom2(i-1:1:i,:) = chrom2(i:-1:i-1,:);
                cchrom2(i-1:1:i,:) = cchrom2(i:-1:i-1,:);
            end
        end
    end
end
%排序重组
TT = [T1;T2];
chrom = [chrom1;chrom2];
cchrom = [cchrom1;cchrom2];
           
function y=yongji(H)
%计算拥挤度
y1=H(:,1);
y2=H(:,2);

[yy1,a1]=sort(y1);
[yy2,a2]=sort(y2);

L=[];
L=[1 1];
for i=2:length(yy1)-1
    L=[L;(yy1(i+1,1)-yy1(i-1,1))/(max(yy1)-min(yy1)),(yy2(i+1,1)-yy2(i-1,1))/(max(yy2)-min(yy2))];
end
L=[L;1 1];
L=[L(a1,1),L(a2,2)];
y=sum(L,2);
end
           

继续阅读