基于matlab-Sheffield工具箱实现遗传算法

发表于 2023-07-08 本文字数 12k 字阅读时长 10 分钟

# 一、Sheffield 工具箱介绍

谢菲尔德（ Sheffield ）遗传算法工具箱是英国谢菲尔德大学开发的遗传算法工具箱，能够为遗传算法研究者和初次使用遗传算法的用户提供了广泛多样的实用函数。但该工具箱不支持多目标优化。

# 二、Sheffield 工具箱的使用方法

（1）下载 Sheffield 工具箱，工具箱已保存至百度网盘进行备份（我的网盘 / 程序 /matlab 遗传算法工具箱 / 篇 1/gatbx）；
（2）将工具箱文件夹（gatbx 文件夹）放在 matlab 的 toolbox 中；
（3）将工具箱添加到 matlab 的搜索路径，添加的代码为 addpath(genpath('C:\Program Files\Polyspace\R2020a\toolbox\gatbx-toolbox')); ，或在 matlab 主页面，依次点击主页→设置路径→添加并包含子文件夹→选择 gabtx 文件夹；
（4）检查是否安装成功，可在命令行窗口输入 ver('gatbx') 进行查看，若显示如下信息，则证明安装成功；

	>> ver('gatbx')
	------------------------------------------------------------------------------------------------
	MATLAB 版本: 9.12.0.1884302 (R2022a)
	MATLAB 许可证编号: 968398
	操作系统: Microsoft Windows 10 家庭中文版 Version 10.0 (Build 19045)
	Java 版本: Java 1.8.0_202-b08 with Oracle Corporation Java HotSpot(TM) 64-Bit Server VM mixed mode
	------------------------------------------------------------------------------------------------
	Genetic Algorithm Toolbox 版本 1.2

# 三、Sheffield 工具箱的使用案例

针对 Sheffield 工具箱的使用，共有两个例子可供参考。第一个是工具箱自带的例程，文件位置为我的网盘 / 程序 /matlab 遗传算法工具箱 / 篇 1 / 遗传算法 /example2.m；第二个为凸模型奇异值分解文章在图像去噪算例的应用，程序如下：

	clear all;
	%% 参数设置
	trun_1=14;trun_2=14;trun_3=14; % 区间方法取的阶次
	zhixindu=1; % 取 3 会漏不少，取 4 基本不会漏，大于 500 找到第一个
	rand_state=5; % 噪声大小，越小噪声越大
	img_gen_par=16300;

	%% 图像的读取及初始化
	original= imread('img_50.jpg'); % 读入图像


	A0 = im2double(original);
	A0_1=A0(:,:,1);
	dim=size(A0_1,1);
	[m,n]=size(A0_1); % 矩阵 A 的大小
	for i=1:m
	for j=1:n
	if A0_1(i,j)==0
	A0_1(i,j)=0.01;
	end
	end
	end
	num_sub_1=rank(A0_1); % 奇异值数量
	[U0_1,s0_1,V0_1]=svd(A0_1); % 矩阵 A 的奇异值分解
	for i=1:m
	if U0_1(1,1)>0
	U0_1=-U0_1;
	end
	end
	for i=1:n
	if U0_1(1,1)>0 % 注意这里不是 V0，还是跟上一个循环一样是 U0
	V0_1=-V0_1;
	end
	end


	A0_2=A0(:,:,2);
	[m,n]=size(A0_2); % 矩阵 A 的大小
	for i=1:m
	for j=1:n
	if A0_2(i,j)==0
	A0_2(i,j)=0.01;
	end
	end
	end
	num_sub_2=rank(A0_2); % 奇异值数量
	[U0_2,s0_2,V0_2]=svd(A0_2); % 矩阵 A 的奇异值分解
	for i=1:m
	if U0_2(1,1)>0
	U0_2=-U0_2;
	end
	end
	for i=1:n
	if U0_2(1,1)>0 % 注意这里不是 V0，还是跟上一个循环一样是 U0
	V0_2=-V0_2;
	end
	end


	A0_3=A0(:,:,3);
	[m,n]=size(A0_3); % 矩阵 A 的大小
	for i=1:m
	for j=1:n
	if A0_3(i,j)==0
	A0_3(i,j)=0.01;
	end
	end
	end
	num_sub_3=rank(A0_3); % 奇异值数量
	[U0_3,s0_3,V0_3]=svd(A0_3); % 矩阵 A 的奇异值分解
	for i=1:m
	if U0_3(1,1)>0
	U0_3=-U0_3;
	end
	end
	for i=1:n
	if U0_3(1,1)>0 % 注意这里不是 V0，还是跟上一个循环一样是 U0
	V0_3=-V0_3;
	end
	end

	%% 不确定性设计
	buquedingzhi=0.01;

	cov_A_1=zeros(m*n, 1);
	for i=1:m
	for j=1:n
	% cov_A_1((i-1)*m+j)=(A0_1(i,j)/rand_state)^2;
	cov_A_1((i-1)*m+j)=buquedingzhi;
	end
	end
	M_1=cov_A_1.^-1;
	DA_1=reshape(cov_A_1,m,n)';

	while true
	A_rand_1=normrnd(A0_1,DA_1.^0.5/zhixindu,m,n);
	sum_temp = 0;
	for i_temp = 1:m
	for j_temp = 1:n
	sum_temp = sum_temp + (A_rand_1(i_temp, j_temp)-A0_1(i_temp, j_temp)) * M_1((i_temp-1)n+j_temp) (A_rand_1(i_temp, j_temp)-A0_1(i_temp, j_temp));
	end
	end
	sum_temp
	if sum_temp<=img_gen_par
	break
	end
	end

	[U_rand_1,s_rand_1,V_rand_1]=svd(A_rand_1);
	RA_1=cov_A_1.^0.5;
	disp('找到第一个不确定性矩阵')


	cov_A_2=zeros(m*n, 1);
	for i=1:m
	for j=1:n
	% cov_A_2((i-1)m+j)=(A0_2(i,j)/rand_state)^2/1210;
	cov_A_2((i-1)*m+j)=buquedingzhi;
	end
	end
	M_2=cov_A_2.^-1;
	DA_2=reshape(cov_A_2,m,n)';

	while true
	A_rand_2=normrnd(A0_2,DA_2.^0.5/zhixindu,m,n);
	sum_temp = 0;
	for i_temp = 1:m
	for j_temp = 1:n
	sum_temp = sum_temp + (A_rand_2(i_temp, j_temp)-A0_2(i_temp, j_temp)) * M_2((i_temp-1)n+j_temp) (A_rand_2(i_temp, j_temp)-A0_2(i_temp, j_temp));
	end
	end
	sum_temp
	if sum_temp<=img_gen_par
	break
	end
	end

	[U_rand_2,s_rand_2,V_rand_2]=svd(A_rand_2);
	RA_2=cov_A_2.^0.5;
	disp('找到第二个不确定性矩阵')



	cov_A_3=zeros(m*n, 1);
	for i=1:m
	for j=1:n
	% cov_A_3((i-1)m+j)=(A0_3(i,j)/rand_state)^2/1210;
	cov_A_3((i-1)*m+j)=buquedingzhi;
	end
	end
	M_3=cov_A_3.^-1;
	DA_3=reshape(cov_A_3,m,n)';

	while true
	A_rand_3=normrnd(A0_3,DA_3.^0.5/zhixindu,m,n);
	sum_temp = 0;
	for i_temp = 1:m
	for j_temp = 1:n
	sum_temp = sum_temp + (A_rand_3(i_temp, j_temp)-A0_3(i_temp, j_temp)) * M_3((i_temp-1)n+j_temp) (A_rand_3(i_temp, j_temp)-A0_3(i_temp, j_temp));
	end
	end
	sum_temp
	if sum_temp<=img_gen_par
	break
	end
	end

	[U_rand_3,s_rand_3,V_rand_3]=svd(A_rand_3);
	RA_3=cov_A_3.^0.5;
	disp('找到第三个不确定性矩阵')



	%% 一阶区间
	% 根据公式 52 计算奇异值对不确定性参数的导数 Hs，矩阵大小为 512512,512，第 mn 个元素表示第 n 个奇异值对第 m 个不确定性参数的导数
	Hs_1=zeros(m*n,num_sub_1);
	for k=1:num_sub_1
	for i=1:m
	for j=1:n
	if i==j
	Hs_1((i-1)n+j,k)=U0_1(j,k)V0_1(i,k);
	end
	end
	end
	end

	% 根据公式 56 计算奇异值的半径 Rs，矩阵大小为 2*1
	Rs_1=zeros(num_sub_1);
	for i=1:num_sub_1
	sum_temp = 0;
	for i_temp = 1:m*n
	sum_temp = sum_temp + Hs_1(i_temp,i)cov_A_1(i_temp)Hs_1(i_temp,i);
	end
	Rs_1(i,i)=sqrt(sum_temp);
	end



	% 根据公式 52 计算奇异值对不确定性参数的导数 Hs，矩阵大小为 512512,512，第 mn 个元素表示第 n 个奇异值对第 m 个不确定性参数的导数
	Hs_2=zeros(m*n,num_sub_2);
	for k=1:num_sub_2
	for i=1:m
	for j=1:n
	if i==j
	Hs_2((i-1)n+j,k)=U0_2(j,k)V0_2(i,k);
	end
	end
	end
	end

	% 根据公式 56 计算奇异值的半径 Rs，矩阵大小为 2*1
	Rs_2=zeros(num_sub_2);
	for i=1:num_sub_2
	sum_temp = 0;
	for i_temp = 1:m*n
	sum_temp = sum_temp + Hs_2(i_temp,i)cov_A_2(i_temp)Hs_2(i_temp,i);
	end
	Rs_2(i,i)=sqrt(sum_temp);
	end





	% 根据公式 52 计算奇异值对不确定性参数的导数 Hs，矩阵大小为 512512,512，第 mn 个元素表示第 n 个奇异值对第 m 个不确定性参数的导数
	Hs_3=zeros(m*n,num_sub_3);
	for k=1:num_sub_3
	for i=1:m
	for j=1:n
	if i==j
	Hs_3((i-1)n+j,k)=U0_3(j,k)V0_3(i,k);
	end
	end
	end
	end

	% 根据公式 56 计算奇异值的半径 Rs，矩阵大小为 2*1
	Rs_3=zeros(num_sub_3);
	for i=1:num_sub_3
	sum_temp = 0;
	for i_temp = 1:m*n
	sum_temp = sum_temp + Hs_3(i_temp,i)cov_A_3(i_temp)Hs_3(i_temp,i);
	end
	Rs_3(i,i)=sqrt(sum_temp);
	end

	% disp ("理论奇异值半径与单次蒙卡奇异值半径之差")
	% diag(Rs_1)-abs(diag(s_rand_1-s0_1));

	%% 计算去噪图像

	% best_s_1=s_rand_1(1:trun_1,1:trun_1);
	%
	% cc_1=U_rand_1(:,1:trun_1)best_s_1V_rand_1(:,1:trun_1)';
	% cc_1(find(cc_1>1))=1;
	% cc_1(find(cc_1<0))=0;
	% gaijin=calculate_SINR(cc_1, A0_1)
	%
	% cc_1=U_rand_1(:,1:trun_1)s_rand_1(1:trun_1,1:trun_1)V_rand_1(:,1:trun_1)';
	% cc_1(find(cc_1>1))=1;
	% cc_1(find(cc_1<0))=0;
	% zaosheng=calculate_SINR(cc_1, A0_1)
	%
	% cc_1=U_rand_1(:,1:trun_1)s0_1(1:trun_1,1:trun_1)V_rand_1(:,1:trun_1)';
	% cc_1(find(cc_1>1))=1;
	% cc_1(find(cc_1<0))=0;
	% yuanshi=calculate_SINR(cc_1, A0_1)
	%
	% s_temp=s_rand_1(1:trun_1,1:trun_1);
	% s_temp(1,1)=s0_1(1,1);
	% cc_1=U_rand_1(:,1:trun_1)s_tempV_rand_1(:,1:trun_1)';
	% cc_1(find(cc_1>1))=1;
	% cc_1(find(cc_1<0))=0;
	% gaibian1=calculate_SINR(cc_1, A0_1)
	%
	% s_temp=s0_1(1:trun_1,1:trun_1);
	% cc_1=U_rand_1(:,1:trun_1)s_tempV_rand_1(:,1:trun_1)';
	% cc_1(find(cc_1>1))=1;
	% cc_1(find(cc_1<0))=0;
	% gaibiansuoyou=calculate_SINR(cc_1, A0_1)


	%% 遗传算法

	NIND=40; % 个体数目
	MAXGEN=50; % 最大遗传代数
	PRECI=20; % 变量的二进制位数
	GGAP=0.95; % 代沟
	px=0.7; % 交叉概率
	pm=0.01; % 变异概率
	trace_1=zeros(trun_1+1,MAXGEN); % 寻优结果的初始值
	Chrom=crtbp(NIND,PRECI*trun_1); % 初始种群
	% FieldD=[PRECI PRECI;lbx lby;ubx uby;1 1;0 0;1 1;1 1]; % 区域描述器
	FieldD=zeros(7,trun_1);
	for i=1:trun_1
	FieldD(1,i)=PRECI;
	FieldD(2,i)=s0_1(i,i)-Rs_1(i,i);
	FieldD(3,i)=s0_1(i,i)+Rs_1(i,i);
	FieldD(4,i)=1;FieldD(5,i)=0;FieldD(6,i)=1;FieldD(7,i)=1;
	end


	gen=0; % 代计数器
	XY=bs2rv(Chrom,FieldD); % 计算初始种群的十进制转换
	ObjV=zeros(size(XY,1),1);
	for i=1:size(XY,1)
	best_s_1=diag(XY(i,:)');
	cc_1=U_rand_1(:,1:trun_1)best_s_1V_rand_1(:,1:trun_1)';
	cc_1(find(cc_1>1))=1;
	cc_1(find(cc_1<0))=0;
	ObjV(i,1)=calculate_SINR(cc_1, A0_1);
	end


	while gen<MAXGEN
	FitnV=ranking(-ObjV); % 分配适应度值
	SelCh=select('sus',Chrom,FitnV,GGAP); % 选择
	SelCh=recombin('xovsp',SelCh,px); % 重组
	SelCh=mut(SelCh,pm); % 变异
	XY=bs2rv(SelCh,FieldD); % 子代个体的十进制转换

	ObjVSel=zeros(size(XY,1),1);
	for i=1:size(XY,1)
	best_s_1=diag(XY(i,:)');
	cc_1=U_rand_1(:,1:trun_1)best_s_1V_rand_1(:,1:trun_1)';
	cc_1(find(cc_1>1))=1;
	cc_1(find(cc_1<0))=0;
	ObjVSel(i,1)=calculate_SINR(cc_1, A0_1);
	end


	[Chrom,ObjV]=reins(Chrom,SelCh,1,1,ObjV,ObjVSel); % 重插入子代到父代，得到新种群
	XY=bs2rv(Chrom,FieldD);
	gen=gen+1; % 代计数器增加
	% 获取每代的最优解及其序号，Y 为最优解，I 为个体的序号
	[Y,I]=max(ObjV);
	trace_1(1:trun_1,gen)=XY(I,:); % 记下每代的最优值
	trace_1(trun_1+1,gen)=Y; % 记下每代的最优值
	end

	optimal_line=trace_1(trun_1+1,:)';






	trace_2=zeros(trun_2+1,MAXGEN); % 寻优结果的初始值
	Chrom=crtbp(NIND,PRECI*trun_2); % 初始种群
	% FieldD=[PRECI PRECI;lbx lby;ubx uby;1 1;0 0;1 1;1 1]; % 区域描述器
	FieldD=zeros(7,trun_2);
	for i=1:trun_2
	FieldD(1,i)=PRECI;
	FieldD(2,i)=s0_2(i,i)-Rs_2(i,i);
	FieldD(3,i)=s0_2(i,i)+Rs_2(i,i);
	FieldD(4,i)=1;FieldD(5,i)=0;FieldD(6,i)=1;FieldD(7,i)=1;
	end


	gen=0; % 代计数器
	XY=bs2rv(Chrom,FieldD); % 计算初始种群的十进制转换
	ObjV=zeros(size(XY,1),1);
	for i=1:size(XY,1)
	best_s_2=diag(XY(i,:)');
	cc_2=U_rand_2(:,1:trun_2)best_s_2V_rand_2(:,1:trun_2)';
	cc_2(find(cc_2>1))=1;
	cc_2(find(cc_2<0))=0;
	ObjV(i,1)=calculate_SINR(cc_2, A0_2);
	end


	while gen<MAXGEN
	FitnV=ranking(-ObjV); % 分配适应度值
	SelCh=select('sus',Chrom,FitnV,GGAP); % 选择
	SelCh=recombin('xovsp',SelCh,px); % 重组
	SelCh=mut(SelCh,pm); % 变异
	XY=bs2rv(SelCh,FieldD); % 子代个体的十进制转换

	ObjVSel=zeros(size(XY,1),1);
	for i=1:size(XY,1)
	best_s_2=diag(XY(i,:)');
	cc_2=U_rand_2(:,1:trun_2)best_s_2V_rand_2(:,1:trun_2)';
	cc_2(find(cc_2>1))=1;
	cc_2(find(cc_2<0))=0;
	ObjVSel(i,1)=calculate_SINR(cc_2, A0_2);
	end


	[Chrom,ObjV]=reins(Chrom,SelCh,1,1,ObjV,ObjVSel); % 重插入子代到父代，得到新种群
	XY=bs2rv(Chrom,FieldD);
	gen=gen+1; % 代计数器增加
	% 获取每代的最优解及其序号，Y 为最优解，I 为个体的序号
	[Y,I]=max(ObjV);
	trace_2(1:trun_2,gen)=XY(I,:); % 记下每代的最优值
	trace_2(trun_2+1,gen)=Y; % 记下每代的最优值
	end

	optimal_line=trace_2(trun_2+1,:)';






	trace_3=zeros(trun_3+1,MAXGEN); % 寻优结果的初始值
	Chrom=crtbp(NIND,PRECI*trun_3); % 初始种群
	% FieldD=[PRECI PRECI;lbx lby;ubx uby;1 1;0 0;1 1;1 1]; % 区域描述器
	FieldD=zeros(7,trun_3);
	for i=1:trun_3
	FieldD(1,i)=PRECI;
	FieldD(2,i)=s0_3(i,i)-Rs_3(i,i);
	FieldD(3,i)=s0_3(i,i)+Rs_3(i,i);
	FieldD(4,i)=1;FieldD(5,i)=0;FieldD(6,i)=1;FieldD(7,i)=1;
	end


	gen=0; % 代计数器
	XY=bs2rv(Chrom,FieldD); % 计算初始种群的十进制转换
	ObjV=zeros(size(XY,1),1);
	for i=1:size(XY,1)
	best_s_3=diag(XY(i,:)');
	cc_3=U_rand_3(:,1:trun_3)best_s_3V_rand_3(:,1:trun_3)';
	cc_3(find(cc_3>1))=1;
	cc_3(find(cc_3<0))=0;
	ObjV(i,1)=calculate_SINR(cc_3, A0_3);
	end


	while gen<MAXGEN
	FitnV=ranking(-ObjV); % 分配适应度值
	SelCh=select('sus',Chrom,FitnV,GGAP); % 选择
	SelCh=recombin('xovsp',SelCh,px); % 重组
	SelCh=mut(SelCh,pm); % 变异
	XY=bs2rv(SelCh,FieldD); % 子代个体的十进制转换

	ObjVSel=zeros(size(XY,1),1);
	for i=1:size(XY,1)
	best_s_3=diag(XY(i,:)');
	cc_3=U_rand_3(:,1:trun_3)best_s_3V_rand_3(:,1:trun_3)';
	cc_3(find(cc_3>1))=1;
	cc_3(find(cc_3<0))=0;
	ObjVSel(i,1)=calculate_SINR(cc_3, A0_3);
	end


	[Chrom,ObjV]=reins(Chrom,SelCh,1,1,ObjV,ObjVSel); % 重插入子代到父代，得到新种群
	XY=bs2rv(Chrom,FieldD);
	gen=gen+1; % 代计数器增加
	% 获取每代的最优解及其序号，Y 为最优解，I 为个体的序号
	[Y,I]=max(ObjV);
	trace_3(1:trun_3,gen)=XY(I,:); % 记下每代的最优值
	trace_3(trun_3+1,gen)=Y; % 记下每代的最优值
	end

	optimal_line=trace_3(trun_3+1,:)';




	%% 绘制对比图像
	best_s_1=diag(trace_1(1:trun_1,MAXGEN));
	cc_1_c=U_rand_1(:,1:trun_1)best_s_1V_rand_1(:,1:trun_1)';
	cc_1_c(find(cc_1_c>1))=1;
	cc_1_c(find(cc_1_c<0))=0;

	cc_1=U_rand_1(:,1:trun_1)s0_1(1:trun_1,1:trun_1)V_rand_1(:,1:trun_1)';
	cc_1(find(cc_1>1))=1;
	cc_1(find(cc_1<0))=0;



	best_s_2=diag(trace_2(1:trun_2,MAXGEN));
	cc_2_c=U_rand_2(:,1:trun_2)best_s_2V_rand_2(:,1:trun_2)';
	cc_2_c(find(cc_2_c>1))=1;
	cc_2_c(find(cc_2_c<0))=0;

	cc_2=U_rand_2(:,1:trun_2)s0_2(1:trun_2,1:trun_2)V_rand_2(:,1:trun_2)';
	cc_2(find(cc_2>1))=1;
	cc_2(find(cc_2<0))=0;




	best_s_3=diag(trace_3(1:trun_3,MAXGEN));
	cc_3_c=U_rand_3(:,1:trun_3)best_s_3V_rand_3(:,1:trun_3)';
	cc_3_c(find(cc_3_c>1))=1;
	cc_3_c(find(cc_3_c<0))=0;

	cc_3=U_rand_3(:,1:trun_3)s0_3(1:trun_3,1:trun_3)V_rand_3(:,1:trun_3)';
	cc_3(find(cc_3>1))=1;
	cc_3(find(cc_3<0))=0;


	figure(1)
	subplot(1,2,1)
	imshow(A0)
	% title (' 原始图像 ')
	subplot(1,2,2)
	A_rand(:,:,1) = A_rand_1;
	A_rand(:,:,2) = A_rand_2;
	A_rand(:,:,3) = A_rand_3;
	imshow(A_rand)
	% title (' 去噪前图像 ')

	figure(2)
	subplot(1,2,1)
	cc(:,:,1) = cc_1;
	cc(:,:,2) = cc_2;
	cc(:,:,3) = cc_3;
	imshow(cc)
	% title (' 常规降噪 ')
	subplot(1,2,2)
	cc_c(:,:,1) = cc_1_c;
	cc_c(:,:,2) = cc_2_c;
	cc_c(:,:,3) = cc_3_c;
	imshow(cc_c)
	% title (' 凸模型降噪 ')

	disp("常规方法的指标值为：")
	disp([calculate_SINR(cc_1, A0_1),calculate_SINR(cc_2, A0_2),calculate_SINR(cc_3, A0_3)])
	disp("凸模型方法的指标值为：")
	disp([calculate_SINR(cc_1_c, A0_1),calculate_SINR(cc_2_c, A0_2),calculate_SINR(cc_3_c, A0_3)])
	disp("性能上升百分比为")
	disp([(calculate_SINR(cc_1_c, A0_1)-calculate_SINR(cc_1, A0_1))/calculate_SINR(cc_1_c, A0_1),...
	(calculate_SINR(cc_2_c, A0_2)-calculate_SINR(cc_2, A0_2))/calculate_SINR(cc_2_c, A0_2),...
	(calculate_SINR(cc_3_c, A0_3)-calculate_SINR(cc_3, A0_3))/calculate_SINR(cc_3_c, A0_3)])
	disp("性能上升平均百分比为")
	disp(((calculate_SINR(cc_1_c, A0_1)-calculate_SINR(cc_1, A0_1))/calculate_SINR(cc_1_c, A0_1)+...
	(calculate_SINR(cc_2_c, A0_2)-calculate_SINR(cc_2, A0_2))/calculate_SINR(cc_2_c, A0_2)+...
	(calculate_SINR(cc_3_c, A0_3)-calculate_SINR(cc_3, A0_3))/calculate_SINR(cc_3_c, A0_3))/3)

	% figure(2)
	% subplot(1,2,1)
	% cc_1_l=U_rand_1(:,1:trun_1)s0_1(1:trun_1,1:trun_1)V_rand_1(:,1:trun_1)';
	% cc_1_l(find(cc_1_l>1))=1;
	% cc_1_l(find(cc_1_l<0))=0;
	% imshow(cc_1_l)
	% title (' 理论最优降噪 ')
	%
	% subplot(1,2,2)
	% imshow(cc_1_c)
	% title (' 凸模型降噪 ')

最优化编程技巧 matlab