遗传算法优化前几代之后开始出错,前几代的最优解看起来都比较正常,主程序单独赋随机值也可以运行出来,不会提示索引出错,是代码不稳定的问题嘛?
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天,您在需要使用的时候【私信】联系我,我会为您补发。