遗传算法求解多目标优化

遗传算法优化前几代之后开始出错,前几代的最优解看起来都比较正常,主程序单独赋随机值也可以运行出来,不会提示索引出错,是代码不稳定的问题嘛?

function [a,b,c,d]=GA()%输入
% 定义变量取值范围
% global C_p C_cm C_pm C_f C_v C_I C_r T_c T_p T_I u0
% C_p=10000;C_cm=1000;C_pm=10000;C_f=2;C_v=4;C_I=977.4;C_r=50;
% T_c=5;T_p=4;T_I=5/60;u0=80;%定义全局变量,输入时间及成本参数,成本单位元,时间单位h
% Field = [tau N kc n];   %各个变量范围
Field =[1 0 1 1;20 20 20 15];   %各个变量范围(约束)
%%定义遗传算法参数
Dim =4;                % 变量维数
PopNum = 10;            % 种群大小(Number of individuals)N一般取20-200
MaxGen = 10;            % 最大遗传代数(Maximum number of generations)一般取值100-1000
Ggap = 0.8;             % 代沟(Generation gap)每代中种群替换的比例
Pc = 0.8;               %交叉概率,不可太偏大或偏小,一般取0.4-0.99
Pm = 0.2;               %变异概率,变异率太高会导致随机性过大,一般取0.005-0.01
SEL_F = 'sus';          % 选择函数名
OBJ_F = 'EC';      % 目标函数名(目标函数)
Pop = Init(Dim,PopNum,Field); % 产生初始种群
PopHis = []; % 记录历史种群中的所有解
ObjHis = []; % 记录历史解的适应值
BestVal = zeros(MaxGen,1);   % 记录每一代的最优适应值
BestSol = zeros(MaxGen,Dim); % 记录每一代的最优解
MeanVal = zeros(MaxGen,1); % 记录每一代的均值
gen = 1; %第一代
[objV,PopHis,ObjHis] = evaluate(OBJ_F,Pop,PopHis,ObjHis); %计算种群的适应值
[BestVal,BestSol] = keepBest(gen,Pop,objV,BestVal,BestSol); %保留第一代最优值
MeanVal(gen) = mean(objV(:,1));
% fprintf('bestVal(%d) = %f  bestSol(%d): T = %f Xo = %d Xp = %d \n',gen,BestVal(gen),gen,BestSol(gen,:)');
while gen < MaxGen
    gen = gen + 1 ;
    FitV = ranking(objV(:,1));                    %分配适应度值(Assign fitness values)
    NewPop = select(SEL_F,Pop,FitV,Ggap);    %选择
    NewPop = crossover(NewPop,Pc);           %交叉
    NewPop = mutate(NewPop,Pm,Field);        %变异
    [objChV,PopHis,ObjHis] = evaluate(OBJ_F,NewPop,PopHis,ObjHis);  %计算子种群的适应值
   [Pop,objV] = reins(Pop,NewPop,1,[1 1.0],objV(:,1),objChV(:,1));%重插入
    [BestVal,BestSol] = keepBest(gen,Pop,objV,BestVal,BestSol); %最优保留
    MeanVal(gen) = mean(objV);
%     fprintf('bestVal(%d) = %f  bestSol(%d): T = %f Xo = %d Xp = %d \n',gen,BestVal(gen),gen,BestSol(gen,:)');
end
 a=BestVal(MaxGen);b=BestSol(MaxGen,:);
% fprintf('bestVal(%d) = %f  bestSol(%d): T = %f Xo = %d Xp = %d \n',MaxGen,BestVal(MaxGen),MaxGen,BestSol(MaxGen,:)');
% figure(1);
% grid on;
% hold on;
% % plot(BestVal);
% gen=1:1:MaxGen;
% plot(gen,BestVal,'-b*');
% hold off;
% saveas(1,'Optimization Result');
% ylabel('平均费用率 ECT元/小时');
% xlabel('遗传代数');
end
function Pop = Init(Dim,PopNum,Field)%初始化
   % tau,N,kc,n,uk
Pop = NaN.*ones(PopNum,Dim);%如何编码
Pop(:,1)=(Field(2,1)-Field(1,1))*rand(PopNum,1)+Field(1,1)*ones(PopNum,1); % tau 实数编码;a+(b-a)*rand的形式
%  Pop(:,1)=randi([Field(1,1) Field(2,1)],PopNum,1);
Pop(:,2)=randi([Field(1,2) Field(2,2)],PopNum,1); % N 整数编码,返回的是整数数组,生成的是POPNUM行一列的数组,是生成的随机
% Pop(:,3)=randi([Field(1,3) Field(2,3)],PopNum,1);
Pop(:,3)=(Field(2,3)-Field(1,3))*rand(PopNum,1)+Field(1,3)*ones(PopNum,1); % kc 实数编码
Pop(:,4)=randi([Field(1,4) Field(2,4)],PopNum,1);%n整数编码
global f
Pop = validate(Pop,f); %种群合法化
end
function [objV,PopHis,ObjHis] = evaluate(OBJ_F,Pop,PopHis,ObjHis)  %依次输出 种群中的所有解,对应的适应值(目标函数值)
[PopNum,~] = size(Pop);   %输出种群的      行数  和  列数  
objV = zeros(PopNum,1);
    for i = 1:PopNum   
        index = isExist(Pop(i,:),PopHis);
        if index~=0          %相等
            objV(i,:) = ObjHis(index,:);   %值的对应
%             fprintf('T = %10f \t Xo = %2f\t Xp = %2f \t\t CR = %f \t ECT = %f is exist.\n',Pop(i,:)',objV(i,:)');
        else
%             fprintf('T = %10f \t Xo = %2f \t Xp = %2f \t\t ',Pop(i,:)');
            [objV(i,1)]= feval(OBJ_F,Pop(i,1),Pop(i,2),Pop(i,3),Pop(i,4));
%             fprintf('CR = %f \t ECT = %f  \t\n',objV(i,:)');
            [nsize,~] = size(PopHis);
            PopHis(nsize + 1,:) = Pop(i,:);
            ObjHis(nsize + 1,:) = objV(i,:);
        end
    end
end
function index = isExist(Pop,PopHis)   % 指示种群  和历史种群   (对应值)是否相等
index = 0;
[nsize,~] = size(PopHis);
for i = 1:nsize
    if isequal(Pop,PopHis(i,:))  %比较两个数组是否相等,维数和值要相等  输出1 
        index = i;
        break; 
    end
end
end
function [BestVal,BestSol] = keepBest(gen,Pop,objV,BestVal,BestSol)
[minVal,minIndex] = min(objV(:,1));
BestVal(gen) = minVal                %保留新的最优适应值
BestSol(gen,:) = Pop(minIndex,:)         %保留新的最优解    
if gen >= 2 && minVal > BestVal(gen-1,1)
    BestVal(gen) = BestVal(gen-1);   %保留历史最优适应值
    BestSol(gen,:) = BestSol(gen-1,:);   %保留历史最优解
end
end
function NewPop = crossover(Pop,Pc)
NewPop = Pop;
for i = 1:2:length(Pop)-1
    if i+1 <= length(Pop)
        if Pc > rand % rand returns a scalar.
        pos = randi([0,3],[1,2]); % 随机的两个个交叉点
        while pos(1,1) == pos(1,2)
            pos = randi([0,3],[1,2]);
        end
        ii = min(pos);
        jj = max(pos);
        alpha = 0.3;
        NewPop(i,ii+1:jj) = alpha * Pop(i+1,ii+1:jj) + (1-alpha) * Pop(i,ii+1:jj); % 算术交叉
        NewPop(i+1,ii+1:jj) = (1-alpha) * Pop(i+1,ii+1:jj) + alpha * Pop(i,ii+1:jj);
        end
    end
end
global f
NewPop = validate(NewPop,f); %种群合法化
end
function Pop = mutate(Pop,Pm,Field)
for i = 1:length(Pop)
    if Pm > rand
        pos = randi(3,1);
        if i==1
            Pop(i,pos) = rand .* (Field(2,pos)-Field(1,pos)) + Field(1,pos);
        else
             Pop(i,pos) = floor(rand .* (Field(2,pos)-Field(1,pos))) + Field(1,pos);      
        end
    end
end
global f
Pop = validate(Pop,f); %种群合理化
end
function Pop = validate(Pop,f)
for i = 1:length(Pop)
    Pop(i,4)=round(Pop(i,4));    %让 o  p   确定为整数 
%     Pop(i,3)=round(Pop(i,3));
%     if Pop(i,3) >= f      %  Xp>=f
%         Pop(i,3) =f-1;   
%     end
%     if Pop(i,2) >= Pop(i,3)      % Xo> = Xp
%         Pop(i,2) = Pop(i,3)-1; 
% %         Pop(i,2) = Pop(i,3); 
%     end
end
end


主程序:
function [ECT] = EC(tau,N,kc,n)
%UNTITLED bino
%   此处显示详细说明
%变量定义
% tau=randperm(20,1) 
% N=randperm(20,1) 
% kc=randperm(20,1)
% n=randperm(15,1)
global lamda v gamma0 gamma1 l Ep p0 p1 uk C_p C_cm C_pm C_f C_v C_I C_r T_c T_p T_I u0 miu xgma
% C_p=10000;C_cm=1000;C_pm=10000;C_f=2;C_v=4;C_I=977.4;C_r=50;
 C_p=1000;C_cm=400;C_pm=300;C_f=5;C_v=0.1;C_I=150;C_r=15;
T_c=0.3;T_p=0.3;T_I=5/60;u0=80;%定义全局变量,输入时间及成本参数,成本单位元,时间单位h
lamda=0.3;v=2;p0=0.01;p1=0.08;uk=20;
miu=80;xgma=8;
gamma0=(normcdf(-kc,0,1));%x-bar是正态,np是二项
gamma1=(normcdf(kc-12*sqrt(n),0,1));
% gamma0=1-(binocdf(kc*n*p0,n,p0)-(binocdf(n*p0,n,p0)));%gamma的过高
% gamma1=binocdf(kc*n*p0,n,1-p1)-(binocdf(n*p0,n,1-p1));
l=randperm(50,1);
Ep=linspace(miu,miu+kc*(xgma/(sqrt(n))),N+1);
Xz=binornd(N,p0,[1,l]);%不合格品数目服从二项分布,还有l如何确定的问题LP
syms t;
g=(lamda*v)*((lamda*t)^(v-1))*exp(-(lamda*t)^v);
PT=[];%给PT数组提前分配内存
PT(1)=1;
% 以下是抽样点概率推导


% 情况1
for k=2:1:l
    if k<=N
        Pt1(k)=0;
    end
   if k>N
    Pt1(k)=(PT(k-N))*(normcdf(Ep(2),miu,xgma)-normcdf(Ep(1),miu,xgma));%二项分布的概率计算有问题
   end
%成功
%情况2
       Pt2(k)=(PT(k-1))*(normcdf(Ep(2),miu,xgma)-normcdf(Ep(1),miu,xgma));
%成功

% % 情况3
    for m=2:1:N-1
      if  k<m+1
        Pt3(k)=0;
      end
      if  k>=m+1
        Pt33(m)=(PT(k-m))*(normcdf(Ep(2),miu,xgma)-normcdf(Ep(1),miu,xgma));
        Pt3(k)=sum(Pt33);
      end
    end
    
%计算概率点   
    if N==1
        PT(k)=1;
    end
    if N==2
        PT(k)=Pt1(k)+Pt2(k);
    end
    if N>2
        PT(k)=Pt1(k)+Pt2(k)+Pt3(k);
    end
end
%计算Wi
for i=2:l
 if N==1
    W(i)=i*tau;
 end
    if N>=2
       W(i)=PT(i)*(i*tau);
    end%有问题
end

%S1,l
for i=2:l
    P1(i)=(int(g,t,[W(i-1),W(i)]))*((1-gamma0)^(i-1))*(1-gamma1);
end
PS1=sum(P1,'all');

%S2(循环正确性)
for i=2:l
    for j=i+1:l
    P2(i,j)=(int(g,t,[W(i-1),W(i)]))*((gamma1)^(j-i))*(1-gamma1)*((1-gamma0)^(i-1));
    end
end
PS2=sum(P2,'all');

%S3
for i=1:l
    P3(i)=((1-int(g,t,[0,W(i)])))*((1-gamma0)^(i-1))*(gamma0);
end
PS3=sum(P3,'all');

%S4
for i=2:l
    P4(i)=(int(g,t,[W(i-1),W(i)]))*((gamma1)^(l-i+1))*((1-gamma0)^(i-1));
end
PS4=sum(P4,'all');

%S5
PS5=1-PS1-PS2-PS3-PS4;
%总概率要差不多等于1
P=PS1+PS2+PS3+PS4+PS5;

%S1
for i=2:1:l
    P1(i)=(int(g,t,[W(i-1),W(i)]))*(gamma0)^i;
    C1(i)=(P1(i))*(C_p+i*(C_f+n*C_v)+C_cm+C_I+(C_r*(sum(Xz(1:i)))));
    T1(i)=(P1(i))*(W(i)+T_c+T_I);
end
EC1=sum(C1,'all');
ET1=sum(T1,'all');
%S2
for i=2:1:l
    for j=i+1:1:l
    P2(i,j)=(int(g,t,[W(i-1),W(i)]))*((1-gamma1)^(j-i))*gamma1;
    A=sum(P2);
    C2(j)=(A(j))*(C_p+j*(C_f+n*C_v)+C_cm+C_I+C_r*(sum(Xz(1:j))));
    T2(j)=(A(j))*(W(j)+T_c+T_I);
    end
end
EC2=sum(C2,'all');
ET2=sum(T2,'all');
%S3
for i=1:1:l
    P3(i)=((1-int(g,t,[0,W(i)])))*((1-gamma0)^i)*gamma0;
    C3(i)=(P3(i))*(C_p+C_pm+i*(C_f+n*C_v)+C_r*(sum(Xz(1:i))));
    T3(i)=(P3(i))*(W(i)+T_I);
end
EC3=sum(C3,'all');
ET3=sum(T3,'all');
%S4
for i=2:1:l
    P4(i)=(int(g,t,[W(i-1),W(i)]))*((1-gamma1)^(l-i+1));
    C4(i)=(P4(i))*(C_p+C_cm+C_I+l*(C_f+n*C_v)+C_r*(sum(Xz(1:l))));
end
EC4=sum(C4,'all');
ET4=(uk*tau+T_c+T_I);
%S5
for i=1:1:l
    PS5=1-(PS1+PS2+PS3+PS4);
    C5(i)=PS5*(C_p+C_pm+C_I+i*(C_f+n*C_v)+C_r*(sum(Xz(1:i))));
end
EC5=sum(C5,'all');
ET5=(uk*tau+T_I);
EC=EC1+EC2+EC3+EC4+EC5;
ET=ET1+ET2+ET3+ET4+ET5;
ECT=vpa(EC/ET);
end

数组索引必须为正整数或逻辑值。

出错 EC (line 36)
Pt1(k)=(PT(k-N))*(normcdf(Ep(2),miu,xgma)-normcdf(Ep(1),miu,xgma));%二项分布的概率计算有问题

出错 GA>evaluate (line 76)
[objV(i,1)]= feval(OBJ_F,Pop(i,1),Pop(i,2),Pop(i,3),Pop(i,4));

出错 GA (line 35)
[objChV,PopHis,ObjHis] = evaluate(OBJ_F,NewPop,PopHis,ObjHis); %计算子种群的适应值
时而报错pt1时而报错pt2,删除了可能出错的索引试运行还是有问题,想知道到底是遗传算法的代码出了问题还是主程序写的问题
顺利优化

你好,我是有问必答小助手,非常抱歉,本次您提出的有问必答问题,技术专家团超时未为您做出解答


本次提问扣除的有问必答次数,将会以问答VIP体验卡(1次有问必答机会、商城购买实体图书享受95折优惠)的形式为您补发到账户。


因为有问必答VIP体验卡有效期仅有1天,您在需要使用的时候【私信】联系我,我会为您补发。