2021數學建模國賽C題思路 生產企業原材料的訂購與運輸 第一版思路 思路開源 已修訂 程式

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)=max(X2')';%供應商最大供應能力，正向指標

Y(:,2)=mean(X2,2);%供應商平均供應能力，正向指標

Y(:,3)=1./std(X2')';%供應商供應的穩定性,對歷史供應資料，做方差，負向指標，可以做個倒數

Y(:,4)=double(X3);%供應商材料經濟效益，正向指標

Y(:,5)=sum(X2,2)./sum(X1,2);%歷史供應量/歷史訂單量比，正向指標

%從企業角度看對供應商的認可情況

Y(:,6)=max(X1')';%供應商最大訂單量，正向指標

Y(:,7)=mean(X1,2);%供應商平均訂單量，正向指標

Y(:,8)=1./std(X1')';%供應商訂單量的穩定性,對歷史供應資料，做方差，負向指標，可以做個倒數

%還有其他指標可自加

Y=log(1+Y);%對指標ln函式映射，以減少資料本身差距，將評分重心轉移至權重

Y=mapminmax(Y',0,1)';

a=[1 5/4 5/4 5/9 5/8 5/6 1 1

4/5 1 1 4/9 4/8 4/6 4/5 4/5

4/5 1 1 4/9 4/8 4/6 4/5 4/5

9/5 9/4 9/4 1 9/8 9/6 9/5 9/5

8/5 2 2 8/9 1 8/6 8/5 8/5

6/5 6/4 6/4 6/9 6/8 1 6/5 6/5

1 5/4 5/4 5/9 5/8 5/6 1 1

1 5/4 5/4 5/9 5/8 5/6 1 1]; %構建判斷矩陣

n1=length(a); %a矩陣的長度

[x,y]=eig(a); %求特征值

lamda=max(diag(y)) %最大特征向量

num=find(diag(y)==lamda);

w0=x(:,num)/sum(x(:,num)) %準則層權向量

ri=[0,0,0.58,0.90,1.12,1.24,1.32,1.41,1.45]; %平均隨機一致性指標

CR1=(lamda-n1)/(n1-1); %一致性指標

CR=CR1/ri(n1); %一致性比例

if CR>=0.10

disp('矩陣沒通過一致性檢驗，請重新調整矩陣');

else

disp('矩陣通過一致性檢驗');

end

f=Y*w0;

[f,paixu]=sort(f,'descend');

[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');

w=[0.0532 0.3138 0.0865 0.0000 0.1076 0.0318 0.1549 0.2523];%通過其他方法求得的權重

ra=tiedrank(Y);

[n,m]=size(ra);%求矩陣維度

RSR=mean(ra,2)/n;

W=repmat(w,[n,1]);%用于復制平鋪矩陣，相當于講w矩陣看成一個元素，形成n×1維的矩陣并用W矩陣記錄

WRSR=sum(ra.*W,2)/n;%加權秩和比

[sWRSR,ind]=sort(WRSR,'descend');%sort為排序函式

p=[1:n]/n;

p(end)=0.99999;

Probit=norminv(p,0,1);

X=[ones(n,1),Probit'];

[ab,abint,r,rint,stats]=regress(sWRSR,X);

WRSRfit=ab(1)+ab(2)*Probit;

f=WRSRfit';

paixu=ind;

[n,m]=size(Y); %n為A矩陣的行數，m為A矩陣的列數

c=sqrt(sum(Y.*Y));

%規范化決策矩陣

d=Y./c;

w=[0.0532 0.3138 0.0865 0.0000 0.1076 0.0318 0.1549 0.2523];%通過其他方法求得的權重

c=w.*d;

cmax=max(c);

cmin=min(c);

for i=1:n

c1=c(i,:)-cmax;

s1(i)=norm(c1);

c2=c(i,:)-cmin;

s2(i)=norm(c2);

T(i)=s2(i)/(s1(i)+s2(i));

end

%排名

[~,pm]=sort(T,'descend');

disp('評分結果，評磁區間[0,1]')

disp(T)

disp('方案排名')

disp(pm)

X=Y;

mapping.mean = mean(X, 1);

X = X - repmat(mapping.mean, [size(X, 1) 1]);%去均值

C = cov(X);%協方差矩陣

C(isnan(C)) = 0;

C(isinf(C)) = 0;

[M, lambda] = eig(C);%求C矩陣特征向量及特征值

[lambda, ind] = sort(diag(lambda), 'descend');%排序

lambda=lambda./sum(lambda);

L=lambda;

lambda=cumsum(lambda);

mapping.lambda = lambda;

k=find(lambda>0.95);

M = M(:,ind(1:k(1)));%%取前k列

mappedX = X * M;%降維后的X

mappedX=mapminmax(mappedX',0,1)';

L=L(1:k(1))./sum(L(1:k(1)));%貢獻率歸一化

f=mappedX*L;

[f,paixu]=sort(f,'descend');

w1=[0.0532 0.3138 0.0865 0.0000 0.1076 0.0318 0.1549 0.2523];%通過其他方法求得的權重

[bestx,f]=youhua(w1,Y);

function [bestx,f1]=youhua(w1,p)

%% 啟發式演算法+投影尋蹤

%可調引數

sub=w1.*0.5; %自變數下限

up=w1.*1.5; %自變數上限

opt=1;%-1為最小化，1為最大化

% 程式為最小化尋優，如果是最大化尋優更改排序部分程式即可

n=length(sub); %自變數個數

num=100*n^2; %種群大小

if num>10000

num=10000;

end

det=10+10*n;%迭代次數

k=1.5+0.1*n;%k為最大環繞圈數

R=0.01*n;%當鯊魚進入獵物該范圍，則直接對獵物位置進行逼近

Mc=(up-sub).*0.1;%獵物行動最大范圍

x=[];

f=[];

for i=1:num

x(i,:)=sub+(up-sub).*rand(1,n); %初始化位置

f(i,1)=Target(p,x);

end

%以最小化為例

if opt==-1

[bestf,a]=min(f);%記錄歷史最優值

bestx=x(a,:);%記錄歷史最優解

elseif opt==1

[bestf,a]=max(f);%記錄歷史最優值

bestx=x(a,:);%記錄歷史最優解

end

trace=[];

trace(1)=bestf;

xx=[];ff=[];

for ii=1:det

%獵物躲避,蒙特卡洛模擬周圍1000次，并選擇最佳的點作為下一逃跑點

d=[];

d=pdist2(bestx,x);

d=sort(d);

z=exp(-d(2)/mean(Mc));%獵物急躁系數

z=max(z,0.1);

dx=[];

yx=[];

bestt=[];

bestc=[];

for i=1:1000

m=[];

for k=1:n

m=[m,(-1)^randi([1,2])];

end

xd=[];

xd=bestx+Mc.*z.*((det-ii)/det).*rand(1,n).*m;

xd=max(sub,xd);

xd=min(up,xd);

dx=[dx;xd];%(det-ii)/det表示隨著追捕，獵物可逃竄的范圍越來越小

yx=[yx;Target(p,xd)];

end

if opt==-1

[bestt,a]=min(yx);%選擇最佳逃跑點

bestc=dx(a,:);

if bestt<bestf

bestf=bestt;

bestx=bestc;

end

elseif opt==1

[bestt,a]=max(yx);%選擇最佳逃跑點

bestc=dx(a,:);

if bestt>bestf

bestf=bestt;

bestx=bestc;

end

end

%鯨魚追捕

for i=1:num

%更新公式

if sqrt(sum((x(i,:)-bestx).^2))<=R

xx(i,:)=x(i,:)+rand.*(x(i,:)-bestx);

xx(i,:)=max(sub,xx(i,:));

xx(i,:)=min(up,xx(i,:));

ff(i,1)=Target(p,xx(i,:));

else

xx(i,:)=x(i,:)+real(shayu(x(i,:)-bestx,k));

xx(i,:)=max(sub,xx(i,:));

xx(i,:)=min(up,xx(i,:));

ff(i,1)=Target(p,xx(i,:));

end

end

%引入上一代進行排序,并重新分配角色

F=[bestf;f;ff];

X=[bestx;x;xx];

if opt==-1

[F,b]=sort(F,'ascend');%按小到大排序

elseif opt==1

[F,b]=sort(F,'descend');%按大到大排序

end

X=X(b,:);

f=F(1:num,:);

x=X(1:num,:);

if opt==-1

[bestf,a]=min(f);%記錄歷史最優值

elseif opt==1

[bestf,a]=max(f);%記錄歷史最優值

end

bestx=x(a,:);%記錄歷史最優解

trace=[trace;bestf];

end

bestx=bestx./sum(bestx);

disp(["投影尋蹤計算權重結果為：",string(bestx)])

f1=[p*bestx'];

function f=shayu(X,K)

%X為鯨魚坐標減去獵物坐標

%k為最大圍繞圈數

n=length(X);

Y=[];

y=[];

costheta=[];

theta=[];

k=[];

for i=1:n-1

Y=[X(1:i),0];

y(i)=sqrt(sum(X(1:i+1).^2));

%計算角度（圈數）

costheta(i)=(X(1:i+1)*Y')/(sqrt(sum(X(1:i+1).^2))*sqrt(sum(Y.^2)));

if isnan(costheta(i))==1

costheta(i)=1;

end

if X(i+1)>=0

theta(i)=acos(costheta(i))/(2*pi);

else

theta(i)=2-acos(costheta(i))/(2*pi);

end

theta(i)=theta(i)*2*pi;

%自適應調節k

if y(i)>=10

k(i)=K*exp(-2);

else

k(i)=K*exp(-y(i)/5);

end

end

%位置更新公式，左包圍或右包圍

f=[];

l=[];

yy=y;

ttheta=theta;

l=k(1).*rand;

yy(1)=yy(1).*exp(-l);

ttheta(1)=ttheta(1)+l*2*pi*(-1).^randi([1,2]);

f=[yy(1).*cos(ttheta(1));yy(1).*sin(ttheta(1))];

if n>2

for j=1:n-2

a=mod(j,2);

if a==1

l=k(j+1).*rand;

yy(j+1)=yy(j+1).*exp(-l);

ttheta(j+1)=ttheta(j+1)+l*2*pi*(-1).^randi([1,2]);

f=[f.*abs(cos(ttheta(j+1)));yy(j+1).*sin(ttheta(j+1))];

elseif a==0

l=k(j+1).*rand;

yy(j+1)=yy(j+1).*exp(-l);

ttheta(j+1)=ttheta(j+1)+l*2*pi*(-1).^randi([1,2]);

f=[f.*abs(sin(ttheta(j+1)));yy(j+1).*cos(ttheta(j+1))];

end

end

end

f=f';

function y=Target(x,a)

%計算目標函式值

m=[];

n=[];

[m,n]=size(x);

z=[];

for i=1:m

s1=0;

for j=1:n

s1=s1+a(j)*x(i,j);

end

z(i)=s1;

end

Sz=[];

Sz=std(z);%方差

R=[];

R=0.1*Sz;

s3=0;

r=[];

t=[];

u=[];

for i=1:m

for j=1:m

r=abs(z(i)-z(j));

t=R-r;

if t>=0

u=1;

else

u=0;

end

s3=s3+t*u;

end

end

Dz=[];

Dz=s3;

y=[];

y=Sz*Dz;

end

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

aa=randperm(402);aa=aa(1:50);

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表格排好格式，直接復制過去即可

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

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];

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

aa=randperm(402);aa=aa(1:50);

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,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