【Matlab】智能优化算法_亨利气体溶解度优化算法HGSO
- 1.背景介绍
- 2.数学模型
- 2.1 亨利定律
- 2.2 HGSO
- 3.文件结构
- 4.伪代码
- 5.详细代码及注释
- 5.1 Create_Groups.m
- 5.2 Evaluate.m
- 5.3 fun_checkpoisions.m
- 5.4 fun_getDefaultOptions.m
- 5.5 HGSO.m
- 5.6 main.m
- 5.7 sumsqu.m
- 5.8 update_positions.m
- 5.9 update_variables.m
- 5.10 worst_agents.m
- 6.运行结果
- 7.参考文献
1.背景介绍
以亨利定律的行为为基础，为HGSO提供了启示。亨利定律，由J.W.亨利于1800年首次提出。总的来说，在特定的温度或压力下，能够溶解在特定量的溶剂中的最大溶质量称为溶解度。因此，HGSO的行为受到了亨利定律的启发。根据等式（1）至（5），亨利定律可用于确定低溶解度气体在液体中的溶解度。
此外，温度和压力是影响溶解度的两个因素；在高温下，固体变得更易溶解，而气体则不易溶解。对于压力，气体的溶解度随着压力的增加而增加。我们的研究涉及气体的溶解度，如图1所示。
2.数学模型
2.1 亨利定律
1803年，威廉·亨利制定了亨利定律，即气体定律。亨利定律规定：在恒定温度下，溶解在给定类型和体积的液体中的给定气体的量与该气体与该液体平衡的分压成正比。因此，亨利定律高度依赖于温度，并表明气体的溶解度（Sg）与气体的分压（Pg）成正比，如以下方程所示：
其中H是亨利常数，该常数在给定温度下对于给定的气体-溶剂组合是特定的，并且气体的分压由Pg表示。
此外，还必须考虑温度对亨利定律常数的影响。亨利定律常数随系统温度的变化而变化，可以用Van’t Hoff方程描述如下：
其中∇sol E是溶解焓，R是气体常数，A和B是H的T依赖性的两个参数。因此，方程（1）可以积分如下：
其中H是参数a和B的函数，其中a和B是H的T依赖性的两个参数。或者，可以在参考温度T=298.15K下基于H创建表达式。
Van’t Hoff方程是有效的，当∇sol E是一个常数时，因此，方程（4）可以重新表述如下：
2.2 HGSO
步骤1：初始化过程。气体的数量（种群大小N）和气体的位置基于以下等式初始化：
其中，种群N中第i个气体的位置由X（i）表示，r是0和1之间的随机数，Xmin、Xmax是问题的边界，t是迭代时间。气体i的数量、类型j的亨利常数值（Hj（t））、簇j中气体i的分压Pi，j和类型j（Ci）的∇solE/R常数值使用以下方程初始化：
其中，l1、l2和l3分别定义为值等于（5E−02、100和1E−02）的常数。
步骤2：聚类。种群因子被划分为等同于气体类型数量的相等簇。每个星团都有相似的气体，因此具有相同的亨利常数值（Hj）。
步骤3：评估。对每个簇j进行评估，以从其类型中的其他簇中识别出达到最高平衡状态的最佳气体。然后，对气体进行排序，以获得整个群体中的最佳气体。
步骤4：更新亨利系数。亨利系数根据以下方程式进行更新：
其中，Hj是簇j的亨利系数，T是温度，Tθ是常数，等于298.15，iter是迭代总数。
步骤5：更新溶解度。溶解度根据以下方程式进行更新：
其中Si，j是气体i在簇j中的溶解度，Pi，j是簇j中气体i的分压，K是常数。
步骤6：更新位置。职位更新如下：
其中气体i在簇j中的位置表示为X（i，j），并且r和t分别是随机常数和迭代时间。X（i，best）是簇j中最好的气体i，而Xbest是群中最好的气。此外，γ是团簇i中气体j与其团簇中气体相互作用的能力，α是其他气体对团簇j中气体i的影响，等于1，β是常数。F（i，j）是簇j中气体i的适应度，相反，Fbest是整个系统中最佳气体的适应度。F是改变搜索代理方向并提供分集=±的标志。
X（i，best）和Xbest是负责平衡勘探和开发能力的两个参数。具体而言，X（i，best）是簇j中最好的气体i，而Xbest是群中最好的气。
步骤7：逃离局部最优。此步骤用于脱离局部最优。使用以下等式对最差代理的数量（Nw）进行排序和选择：
其中N是搜索代理的数量。
步骤8：更新最差代理的位置。
其中，G（i，j）是气体i在簇j中的位置，r是随机数，Gmin，Gmax是问题的边界。
最后，给出了该算法的伪代码，包括初始种群大小、种群评估和更新参数。
3.文件结构
Create_Groups.m % 创建组
Evaluate.m % 评估位置
fun_checkpoisions.m % 检查位置
fun_getDefaultOptions.m % 获得默认选项
HGSO.m % 亨利气体溶解度优化算法
main.m % 主函数
sumsqu.m % 归属权
update_positions.m % 更新位置
update_variables.m % 更新变量
worst_agents.m % 获取最差代理
4.伪代码
5.详细代码及注释
5.1 Create_Groups.m
function[Group]=Create_Groups(var_n_gases,var_n_types,X)
N=var_n_gases/var_n_types;
i=1;
for j=1:var_n_types
Group{j}.Position=X(i:i+N,:);
i=j*N+1;
if i+N>var_n_gases
i= j*N;
end
end
end
5.2 Evaluate.m
function [X,best_fit,best_pos] = Evaluate(objfunc,var_n_types,var_n_gases, X,Xnew,init_flag)
if init_flag==1
for j=1:var_n_gases/var_n_types
X.fitness(j) = objfunc(X.Position(j,:));
end
else
for j=1:var_n_gases/var_n_types
temp_fit = objfunc(Xnew.Position(j,:));
if temp_fit<X.fitness(j)
X.fitness(j)=temp_fit;
X.Position(j,:)= Xnew.Position(j,:);
end
end
end
[best_fit,index_best]=min(X.fitness(:));
best_pos=X.Position(index_best,:);
end
5.3 fun_checkpoisions.m
function Group=fun_checkpoisions(dim,Group,var_n_gases,var_n_types,var_down,var_up)
Lb=var_down*ones(1,dim);
% Upper bounds
Ub=var_up*ones(1,dim);
for j=1:var_n_types
for i=1:var_n_gases/var_n_types
isBelow1 = Group{j}.Position(i,:) < Lb;
isAboveMax = (Group{j}.Position(i,:) > Ub);
if isBelow1 == true
Group{j}.Position(i,:) =Lb;
elseif find(isAboveMax== true)
Group{j}.Position(i,:) = Ub;
end
end
end
end
5.4 fun_getDefaultOptions.m
function [var_n_gases,var_n_types,var_niter] = fun_getDefaultOptions()
var_n_gases = 35; % The swarm size.
var_n_types = 5; % The number of group.
var_niter = 1000; % The number of iterations.
end
5.5 HGSO.m
function [vec_Xbest, var_Gbest,vec_Gbest_iter] = HGSO(objfunc, dim,var_down,var_up,var_niter,var_n_gases,var_n_types)
if nargin<5
[var_n_gases,var_n_types,var_niter]=fun_getDefaultOptions();
end
%constants in eq (7)
l1=5E-03;
l2=100;
l3=1E-02;
%constants in eq (10)
alpha=1;
beta=1;
%constant in eq (11)
M1=0.1;
M2=0.2;
%paramters setting in eq. (7)
K=l1*rand(var_n_types,1);
P=l2*rand(var_n_gases,1);
C=l3*rand(var_n_types,1);
%randomly initializes the position of agents in the search space
X=var_down+rand(var_n_gases,dim)*(var_up-var_down);
%The population agents are divided into equal clusters with the same Henry痴 constant value
Group=Create_Groups(var_n_gases,var_n_types,X);
% Compute cost of each agent
for i = 1:var_n_types
[Group{i},best_fit(i), best_pos{i}] = Evaluate(objfunc,var_n_types,var_n_gases, Group{i},0,1);
end
[var_Gbest, var_gbest] = min(best_fit);
vec_Xbest = best_pos{var_gbest};
for var_iter = 1:var_niter
[S]=update_variables(var_iter,var_niter,K,P,C,var_n_types,var_n_gases);
Groupnew=update_positions(Group,best_pos,vec_Xbest,S,var_n_gases,var_n_types,var_Gbest,alpha,beta,dim);
Groupnew=fun_checkpoisions(dim,Groupnew,var_n_gases,var_n_types,var_down,var_up);
for i = 1:var_n_types
[Group{i},best_fit(i), best_pos{i}] = Evaluate(objfunc,var_n_types,var_n_gases, Group{i},Groupnew{i},0);
Group{i}=worst_agents(Group{i},M1,M2,dim,var_up,var_down,var_n_gases,var_n_types);
end
[var_Ybest, var_index] = min(best_fit);
vec_Gbest_iter(var_iter)=var_Ybest;
if var_Ybest<var_Gbest
var_Gbest=var_Ybest;
vec_Xbest = best_pos{var_index};
end
end
end
5.6 main.m
clc;clear all;
close all;
fitfun = @sumsqu;
dim=30;
var_niter=1000;
Lb=-10;
Ub=10;
% var_n_gases=35;
% var_n_types=5;
[xf,fval,vec_Gbest_iter]=HGSO(fitfun,dim,Lb,Ub);
figure,
semilogy(vec_Gbest_iter,'r')
xlim([0 var_niter]);
5.7 sumsqu.m
function [y] = sumsqu(xx)
d = length(xx);
sum = 0;
for ii = 1:d
xi = xx(ii);
sum = sum + ii*xi^2;
end
y = sum;
end
5.8 update_positions.m
function [Group]=update_positions(Group,best_pos,vec_Xbest,S,var_n_gases,var_n_types,var_Gbest,alpha,beta,var_nvars)
vec_flag=[1,-1];
for i=1:var_n_types
for j=1:var_n_gases/var_n_types;
gama=beta*exp(-(var_Gbest+.05)/(Group{i}.fitness(j)+.05));
flag_index = floor(2*rand()+1);
var_flag=vec_flag(flag_index);
for k=1:var_nvars
Group{i}.Position(j,k)= Group{i}.Position(j,k)+var_flag*rand*gama*(best_pos{i}(k)-Group{i}.Position(j,k))+rand*alpha*var_flag*(S(i)*vec_Xbest(k)-Group{i}.Position(j,k));
end
end
end
end
5.9 update_variables.m
function [S]=update_variables(var_iter,var_niter,K,P,C,var_n_types,var_n_gases)
T=exp(-var_iter/var_niter);
T0= 298.15;
i=1;
N=var_n_gases/var_n_types;
for j=1:var_n_types
K(j)=K(j)*exp(-C(j)*(1/T-1/T0));
S(i:i+N,:)=P(i:i+N,:)*K(j);
i=j*N+1;
if i+N>var_n_gases
i= j*N;
end
end
end
5.10 worst_agents.m
function [X]=worst_agents(X,M1,M2,dim,G_max,G_min,var_n_gases,var_n_types)
%Rank and select number of worst agents eq.(11)
[X_sort,X_index]=sort(X.fitness(:),'descend');
M1N =M1*var_n_gases/var_n_types;
M2N = M2*var_n_gases/var_n_types;
Nw = round((M2N-M1N).*rand(1,1) + M1N);
for k=1:Nw
X.Position(X_index(k),:)=G_min+rand(1,dim)*(G_max-G_min);
end
end
6.运行结果
7.参考文献
[1]Hashim A F,Houssein H E,Mabrouk S M, et al. Henry gas solubility optimization: A novel physics-based algorithm[J]. Future Generation Computer Systems,2019,101©.
