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