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

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

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

---

前言

本文以第二问为例介绍代码：
问题2. 附件3给出了某次任务的飞行器计划耗油速度数据，与飞行器在飞行器坐标系下的理想质心位置数据。根据任务需求，在飞行器始终保持平飞(俯仰角为0)的任务规划过程中，请为飞行器制定该次任务满足条件(1)∼(6)的6个油箱供油策略，使得飞行器每一时刻的质心位置与理想质心位置的欧氏距离的最大值达到最小, i.e,

请给出飞行器飞行过程中6个油箱各自的供油速度曲线和4个主油箱的总供油速度曲线(时间间隔为1s)、以及飞行器瞬时质心与理想质心距离的最大值和4个主油箱的总供油量，并将6个油箱的供油速度数据按时间（每秒一组）先后顺序存入附件6结果表“第二问结果”中。

---

提示：本文以第二问为例介绍代码：

一、主程序

示例：利用遗传算法求解寻优

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_all=xlsread('..\data\附件3-问题2数据.xlsx','发动机耗油速度');
vc=vc_all(:,2);
%读入飞行器理想质心
ci_all=xlsread('..\data\附件3-问题2数据.xlsx','飞行器理想质心数据');
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
%% 开始生成种群
flag=0;flagnumber=0;alpha_vcp=1.1;line_vcp=3;b=ones(1,6);rn=1;addret=0;addv=0;
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(2,:)=[bound(2,2)*b(2)+addret,Vt_left(2),(vcp*alpha_vcp+line_vcp)/IOnumber+addret];
    boundnumber(3,:)=[bound(3,2)*b(3)+addret,Vt_left(3),(vcp*alpha_vcp+line_vcp)/IOnumber+addret];
    boundnumber(4,:)=[bound(4,2)*b(4)+addret,Vt_left(4),(vcp*alpha_vcp+line_vcp)/IOnumber+addret];
    boundnumber(5,:)=[bound(5,2)*b(5)+addret,Vt_left(5),(vcp*alpha_vcp+line_vcp)/IOnumber+addret];
    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(2:5)=ret(2:5)+addv;
    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阶段惩罚');
        addv=0.3;
        line_vcp=line_vcp+1;
        addret=addret+0.5;
        b=[1 1.1 1.1 1.1 1.1 1];
        rn=1.5*rn;
    elseif flagnumber>=110 && flagnumber<120
        %disp('正在寻找初始化种群,并启动2阶段惩罚');
        addv=0.4;
        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阶段惩罚');
        addv=0.5;
        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
        addv=0.8;
        addret=addret+0.5;
        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;

总结

1.运行结果展示：

2.写在最后

运行结果还有不完善的地方请批评指正,若需求其他程序代码请私信我~~~~

感谢观看~~~