# “华为杯”第十七届中国研究生数学建模竞赛F题代码分享

F题：飞行器质心平衡供油策略优化研究

# 一、主程序

clc
clear
%% Step1：parameters
maxgen = 10;%最大代数
sizepop = 100;%种群规模>30
pcross = [0.60];%交叉概率
pmutation = [0.1];%变异概率
lenchrom = [1 1 1 1 1 1];%染色体的长度
bound = [0 1.1;0 1.8;0 1.7;0 1.5;0 1.6;0 1.1];%速度上限
Vt_init = [0.3 1.5 2.1 1.9 2.6 0.8]*850;%初始油量,kg
v_contin =[0 0 0 0 0 0];%发动机供油口开关计数
%读入发动机耗油速度
vc=vc_all(:,2);
%读入飞行器理想质心
ci=ci_all(:,2:4);
vc_6=[0 0 0 0 0 0];%初始耗油速度
io=zeros(1,6);
dnum = 60;
%% 全程计算
for p=1:length(vc)
%更新剩余油量
Vt_init(p,1)=Vt_init(p,1)-vc_6(p,1);
Vt_init(p,2)=Vt_init(p,2)-vc_6(p,2)+vc_6(p,1);
Vt_init(p,3)=Vt_init(p,3)-vc_6(p,3);
Vt_init(p,4)=Vt_init(p,4)-vc_6(p,4);
Vt_init(p,5)=Vt_init(p,5)-vc_6(p,5)+vc_6(p,6);
Vt_init(p,6)=Vt_init(p,6)-vc_6(p,6);
%更新开关量
%% Step2: individual initialization
individuals = struct('fitness',zeros(1,sizepop),'chrom',[]);
avgfitness = [];
bestfitness = [];
bestchrom = [];
%% Step3： initial population
for i=1:sizepop
[individuals.chrom(i,:),dnum]=CodeNew(lenchrom,bound,v_contin,vc(p),Vt_init(p,:),dnum);
x=individuals.chrom(i,:);
individuals.fitness(i)=fun(Vt_init(p,:)/850,x(1),x(2),x(3),x(4),x(5),x(6),ci(p,:));
end
%disp('初始化完成');
%% Step4： Find the best chromosome
[bestfitness,bestindex] = min(individuals.fitness);
bestchrom = individuals.chrom(bestindex,:);
%avgfitness = sum(individuals.fitness)/sizepop;
%% Step5： Record the best fitness and average fitness
trace = [];
%% Step6： Evolution beginning
for i=1:maxgen
individuals = select(individuals,sizepop);%根据适应度选择
avgfitness = sum(individuals.fitness)/sizepop;%平均适应度

individuals.chrom = Cross(pcross,lenchrom,individuals.chrom,sizepop,bound,vc(p),Vt_init(p,:));
individuals.chrom = Mutation(pcross,lenchrom,individuals.chrom,sizepop,bound,[i maxgen],vc(p),Vt_init(p,:));
%     if mod(i,10) == 0
%         individuals.chrom =nonlinear(individuals.chrom,sizepop);
%     end
for j = 1:sizepop
x = individuals.chrom(j,:);
individuals.fitness(j) = fun(Vt_init(p,:)/850,x(1),x(2),x(3),x(4),x(5),x(6),ci(p,:));
end
[newbestfitness,newbestindex] =min(individuals.fitness);
if bestfitness>newbestfitness
bestfitness = newbestfitness;
bestchrom = individuals.chrom(newbestindex,:);
end
avgfitness = sum(individuals.fitness)/sizepop;
trace = [trace;avgfitness bestfitness];
end
%     figure(1);
%     plot(1:maxgen,trace(:,1),'r');
%     hold on;
%     plot(1:maxgen,trace(:,2),'b');
%     hold off;
fprintf(['第',num2str(p),'秒的optimization :\nv1 v2 v3 v4 v5 v6=    ',num2str(bestchrom),'\nminf=    ',num2str(bestfitness),'\n']);
%% 判断计数是否清零
if p==1
[~,sw]=find(bestchrom>0);
v_contin(sw)=v_contin(sw)+1;
else
[~,sw]=find(vc_6(p,:)>0);
for i=1:length(sw)
if bestchrom(sw(i))>0
v_contin(sw(i))=v_contin(sw(i))+1;
elseif bestchrom(sw(i)) == 0
v_contin(sw(i))=0;
end
end
[~,sw]=find(bestchrom>0);
for i=1:length(sw)
if vc_6(p,sw(i))==0
v_contin(sw(i))=1;
end
end
end

%% 更新存储迭代量
Vt_init=[Vt_init;Vt_init(p,:)];
vc_6=[vc_6;bestchrom];
if dnum >60
dnum =dnum-1;
end
end


# 二、子程序代码

## 1.编码代码

function [ret,dnum]=Code(lenchrom,bound,vkey,vcp,Vt_left,dnum)
%本函数将变量编码成染色体，用于随机初始化一个种群
% lenchrom   input : 染色体长度
% bound      input : 变量的取值范围
% ret        output: 染色体的编码值
%% 检验油口打开时间，超过60s关闭后重新规划
[~,sw]=find(vkey>0);
for i=1:length(sw)
if vkey(sw(i))>=60
vkey(sw(i))=0;
end
end
% if (dum<62 && dum>60) || (dum<120)
%     if vkey(1)==0 && vkey(2)~=0
%          vkey(1)=0;
%     elseif vkey(2)==0 && vkey(1)~=0
%          vkey(2)=0;
%     elseif vkey(2)==0 && vkey(1)==0
%         vkey(1)=0;
% end
%% 开始生成种群
while flag==0 || flag==-1 || flag==-2 || flag==-3 || flag==-4
%% 确认目前必须打开的油口
I0 = [0 0 0 0 0 0]';
[a,sw]=find(vkey>0);
I0(sw)=1;
%% 重新规划油口
flag_sw = length(sw);%已经开的口的数量
if sum(I0(2:5)) ==0 && (I0(1)+I0(6))==0
flag_opens25 = randi(2);
if flag_opens25 == 1
[~,posmax]=sort(Vt_left(2:5));
I0(posmax(end)+1)=1;
%             Idex_list = randperm(4)+1;
%             Idex = Idex_list(1);
%             I0(Idex)=1;
elseif flag_opens25 == 2
[~,posmax]=sort(Vt_left(2:5));
I0(posmax(end-1:end)+1)=1;
%             Idex_list = randperm(4)+1;
%             Idex = Idex_list(1:2);
%             I0(Idex)=1;
end
flag_opens16 = randi(2);
if flag_opens16 == 1
I0(1)=1;
elseif flag_opens16 == 2
I0(6)=1;
end
elseif sum(I0(2:5)) ==0 && (I0(1)+I0(6))==1
flag_opens25 = randi(2);
if flag_opens25 == 1
[~,posmax]=sort(Vt_left(2:5));
I0(posmax(end)+1)=1;
%             Idex_list = randperm(4)+1;
%             Idex = Idex_list(1);
%             I0(Idex)=1;
elseif flag_opens25 == 2
[~,posmax]=sort(Vt_left(2:5));
I0(posmax(end-1:end)+1)=1;
%             Idex_list = randperm(4)+1;
%             Idex = Idex_list(1:2);
%             I0(Idex)=1;
end
elseif sum(I0(2:5)) ==1 && (I0(1)+I0(6))==0
[~,posmax]=sort(Vt_left(2:5));
%[~,pos]=find(I0(2:5)==1);
hg = I0(2:5);
if hg(posmax(end)) == 1
I0(posmax(end-1)+1)=1;
elseif hg(posmax(end)) == 0
I0(posmax(end)+1)=1;
end
%         while flag_pos==0
%             Idex_list = randperm(4)+1;
%             Idex = Idex_list(1);
%             if hg(pos) == Idex
%                 flag_pos = 0;
%             else
%                 flag_pos = 1;
%             end
%         end
%         I0(Idex)=1;
flag_opens16 = randi(2);
if flag_opens16 == 1
I0(1)=1;
elseif flag_opens16 == 2
I0(6)=1;
end
elseif sum(I0(2:5)) ==1 && (I0(1)+I0(6))==1
[~,posmax]=sort(Vt_left(2:5));
%[~,pos]=find(I0(2:5)==1);
hg = I0(2:5);
if hg(posmax(end)) == 1
I0(posmax(end-1)+1)=1;
elseif hg(posmax(end)) == 0
I0(posmax(end)+1)=1;
end
%         [~,pos]=find(I0(2:5)==1);
%         flag_pos=0;
%         hg = I0(2:5);
%         while flag_pos==0
%             Idex_list = randperm(4)+1;
%             Idex = Idex_list(1);
%             if hg(pos) == Idex
%                 flag_pos = 0;
%             else
%                 flag_pos = 1;
%             end
%         end
%         I0(Idex)=1;
elseif sum(I0(2:5)) ==2 && (I0(1)+I0(6))==0
flag_opens16 = randi(2);
if flag_opens16 == 1
I0(1)=1;
elseif flag_opens16 == 2
I0(6)=1;
end
end
%% 生成染色体
IOnumber = sum(I0(2:5))+1;
pick=1-rand(1,length(lenchrom));
while min(pick) ==0
pick=rand(1,length(lenchrom));
disp('1');
end
retbound=zeros(6,2);
retbound(:,1)=bound(:,1);
retbound(2:5,1)=vcp/2.1;
%判断油箱剩余油量从而永久关闭油箱
%     if Vt_left(1)<0.00001
%         io(1)=1;%1代表抑制，0代表激活
%         %     elseif Vt_left(3)<0.001
%         %         I0(3)=0;
%         %     elseif Vt_left(4)<0.001
%         %         I0(4)=0;
%     elseif Vt_left(6)<0.00001
%         io(6)=1;
%     end
%vcp宽裕系数
boundnumber=zeros(6,3);
boundnumber(1,:)=[bound(1,2)*0.1,1645.6-Vt_left(2),Vt_left(1)];
boundnumber(6,:)=[bound(6,2)*0.1,2448-Vt_left(5),Vt_left(6)];
retbound(1,2)=min(boundnumber(1,:));
%     retbound(2,2)=min(min(bound(2,2),Vt_left(2)),(vcp*alpha_vcp+line_vcp)/IOnumber);
%     retbound(3,2)=min(min(bound(3,2),Vt_left(3)),(vcp*alpha_vcp+line_vcp)/IOnumber);
%     retbound(4,2)=min(min(bound(4,2),Vt_left(4)),(vcp*alpha_vcp+line_vcp)/IOnumber);
%     retbound(5,2)=min(min(bound(5,2),Vt_left(5)),(vcp*alpha_vcp+line_vcp)/IOnumber);
retbound(2,2)=abs(min(boundnumber(6,:)));
retbound(3,2)=abs(min(boundnumber(6,:)));
retbound(4,2)=abs(min(boundnumber(6,:)));
retbound(5,2)=abs(min(boundnumber(6,:)));
retbound(6,2)=abs(min(boundnumber(6,:)));
ret=retbound(:,1)'+(retbound(:,2)-retbound(:,1))'.*pick; %线性插值，编码结果以实数向量存入ret中
ret=ret.*I0';
flag=test(lenchrom,bound,ret,vcp,Vt_left);     %检验染色体的可行性
if flag ==-1
ret;
end
%I0'
%ret
%flag
if flag==0 || flag==-1 || flag==-2 || flag==-3 || flag==-4
flagnumber=flagnumber+1;
end
%惩罚措施
if 110>flagnumber && flagnumber>=100
%disp('正在寻找初始化种群,并启动1阶段惩罚');
line_vcp=line_vcp+1;
b=[1 1.1 1.1 1.1 1.1 1];
rn=1.5*rn;
elseif flagnumber>=110 && flagnumber<120
%disp('正在寻找初始化种群,并启动2阶段惩罚');
line_vcp=line_vcp+1;
b=[1 1.3 1.3 1.3 1.3 1];
rn=2;
elseif flagnumber>=120 && flagnumber<130
%disp('正在寻找初始化种群,并启动3阶段惩罚');
alpha_vcp=alpha_vcp+0.05;
b=[1 1.5 1.5 1.5 1.5 1];
rn=2;
elseif flagnumber>=130 && flagnumber<200
%disp('正在寻找初始化种群,并启动4阶段惩罚');
alpha_vcp=alpha_vcp+0.1;
b=[1 1.8 1.8 1.8 1.8 1];
rn=2;
elseif flagnumber>=200 && flagnumber<210
disp('正在寻找初始化种群,并启动5阶段惩罚');
line_vcp=line_vcp+2;
rn=3;
elseif flagnumber>=210
disp('正在努力初始化');
flag
b=[1 1.3 1.3 1.3 1.3 1];
line_vcp=line_vcp+1;
rn=1.1*rn;
end

end
%disp('编码完成');


## 2.判断代码

function flag=test(lenchrom,bound,code,vcp,Vt_left)
% lenchrom   input : 染色体长度
% bound      input : 变量的取值范围
% code       output: 染色体的编码值
%x=code; %先解码
flag=1;
%% 判断油口正确性
if flag == 1
%提取油箱供油口
hg=code(:,2:5);
[~,m] =find(code>0);%总供油口数
[~,p] =find(hg>0);%发动机供油口数
if length(m)>3 || length(p) >2
flag=-1;
%disp('油口错误');
%else
%disp('油口正确');
end
end

%% 判断供需条件
if flag == 1
Vx = vcp;
if sum(hg) >Vx
flag=1;
%disp('供需正确');
else
%flag=0;
flag=-2;
end
end

%% 判断速度范围条件
[n,~]=size(code);
if flag ==1
for i=1:n
if code(i)<bound(i,1) || code(i)>bound(i,2)
flag=-3;
disp('速度错误');
end
end
end
%% 判断油箱剩余油量是否为负
[~,m]=size(Vt_left);
if flag ==1
for i=1:m
if Vt_left(i)<0
flag=-4
disp('剩余油量错误');
end
end
end
% if (x(1)<0)&&(x(2)<0)&&(x(3)<0)&&(x(1)>bound(1,2))&&(x(2)>bound(2,2))&&(x(3)>bound(3,2))
%     flag=0;
% end


## 3.交叉代码

function ret = Cross(pcross,lechrom,chrom,sizepop,bound,vcp,Vt_left)
for i=1:sizepop
flag = 0;
%% 寻找位数相同的两个染色体进行交叉
flag_pos=0;flag_posnum=0;
while flag_pos ==0
%先选择交叉的位置
Idex_list = randperm(6);
pos = Idex_list(1);
%筛选出这个位置油口打开的染色体
Ic=[];
for i=1:sizepop
if chrom(i,pos)>0
Ic=[Ic;i];
end
end
%判断个数是否大于2
if length(Ic)<2
flag_pos =0;
else
flag_pos=1;
index=Ic(1:2);
end
%在这些染色体中任意选择前两个进行交叉
flag_posnum=flag_posnum+1;
if flag_posnum>30
disp('无交叉对象');
flag =-1;
break;
end
%disp('3');
end
%     pick =rand(1,2);picknum=0;
%     while prod(pick) == 0%对pick1行2列中的两个数进行相乘，判断是否为0
%         pick =rand(1,2);
%         picknum=picknum+1;
%         if picknum >30
%             disp('Cross中的pick=rand(1,2)寻找失败！');
%             continue;
%         end
%     end
%     index = ceil(pick.*sizepop);
%%是否进行交叉动作
pick = rand;
while pick == 0
pick = rand;
disp('5');
end
if pick>pcross
continue;
end
flag = 0;flagnum=0;
while flag ==0
%         while flag_pos ==0
%             pick = rand;
%             while pick ==0
%                 pick = rand;
%             end
%             pos = ceil(pick.*sum(lechrom));%向无穷大的方向取整
%             if chrom(index(1),pos)>0 && chrom(index(2),pos)>0
%                 flag_pos=1;
%             end
%             flag_posnum=flag_posnum+1;
%             if flag_posnum>4
%                 disp('交叉失败,跳出循环！');
%                 continue;
%             end
%         end
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(lechrom,bound,chrom(index(1),:),vcp,Vt_left);
%flag1
flag2 = test(lechrom,bound,chrom(index(2),:),vcp,Vt_left);
%flag2
if flag1*flag2 == 0
flag = 0;
else
flag =1;
end
flagnum=flagnum+1;
if flagnum>1
disp('交叉环节失败！');
break;
end
end
end
ret = chrom;
end


## 4.变异代码

function ret = Mutation(pmutation,lechrom,chrom,sizepop,bound,pop,vcp,Vt_left)
for i=1:sizepop
flag_pos=0;flag_posnum=0;
while flag_pos ==0
%先选择变异的位置
Idex_list = randperm(6);
pos = Idex_list(1);
%筛选出这个位置油口打开的染色体
Ic=[];
for i=1:sizepop
if chrom(i,pos)>0
Ic=[Ic;i];
end
end
%判断个数是否大于2
if length(Ic)<1
flag_pos =0;
else
flag_pos=1;
end
%在这些染色体中任意选择前两个进行交叉
flag_posnum=flag_posnum+1;
if flag_posnum>30
disp('无变异对象');
flag =-1;
break;
end
%disp('3');
end
pick =rand;
if pick>pmutation
continue;
end
flag =0;flag_number=0;
while flag ==0
pick =rand;
while pick ==0
pick =rand;
end
for j=1;length(Ic);
v=chrom(j,pos);
v1=v-bound(pos,1);
v2=bound(pos,2)-v;
pick=rand;
if pick >0.5
delta = v2*(1-pick^((1-pop(1)/pop(2))^2));
chrom(i,pos)=v+delta;
else
delta = v1*(1-pick^((1-pop(1)/pop(2))^2));
chrom(i,pos) = v-delta;
end
flag=test(lechrom,bound,chrom(i,:),vcp,Vt_left);
flag_number=flag_number+1;
if flag_number>1
disp('变异环节失败！');
break;
end
end
end
end
ret = chrom;
end


## 5.选择代码

function ret = select(individuals,sizepop)
individuals.fitness = 1./(individuals.fitness);
sumfitness = sum(individuals.fitness);
sumf = individuals.fitness./sumfitness;%归一化
index=[];
for i=1:sizepop
pick = rand;picknum=0;
while pick ==0
pick = rand;
picknum=picknum+1;
if picknum >30
disp('Select中的pick寻找失败！');
end
end
for j=1:sizepop
pick = pick-sumf(j);
if pick<0
index = [index j];
break;
end
end
end
individuals.chrom = individuals.chrom(index,:);
individuals.fitness = individuals.fitness(index);
ret = individuals;


# 总结

## 2.写在最后

• 点赞 3
• 评论
• 分享
x

海报分享

扫一扫，分享海报

• 收藏 1
• 手机看

分享到微信朋友圈

x

扫一扫，手机阅读

• 打赏

打赏

Python_Cecil

你的鼓励将是我创作的最大动力

C币 余额
2C币 4C币 6C币 10C币 20C币 50C币
• 一键三连

点赞Mark关注该博主, 随时了解TA的最新博文
10-24

11-11 2603
09-21 469
09-23 3756