MATLAB 智能算法入门学习
matlab常用命令
https://mp.weixin.qq.com/s/_V-jvcUI_CfIHpVZOeogjQ
NO.1 画图类
----------------------------------
plot 最常用的画图命令，画线的时候坐标要注意
subplot 生成子图
plot3 生成线条组成的三维图形
mesh 生成网格状立体图形
meshgrid 生成mesh需要的网格数据点
ezplot 主要用来画隐函数的图像
scatter 主要画散点图
figure 生成一张图
cla 清除轴线图形
title 设置标题
axis 设置坐标范围，也有manual、equal等操作
grid on 生成网格
set 设置，和gca，gcf一起用比较多
xlabel，ylabel 设置坐标轴的名字大小字体
text 在图上进行标注，可以指定位置、颜色等
hold on 保持，可以用来显示多个图形
legend 图例
pause 暂停命令，我常用来做动态效果
xlim 设置坐标范围
gca 返回当前axes 对象的句柄值，经常和set一起用
gcf 返回当前对象的句柄值，经常和set一起用
aexs 创建一个坐标系，GUI中常用
NO.2 运算操作
----------------------------------
sin() 正弦函数
asin 反正弦函数
cos 余弦函数
acos 反余弦函数
tan 正切函数
atan 反正切函数
exp 指数函数
sqrt 开根号
log 对数
abs 取绝对值
sum 向量元素的求和
pi 3.14159那个，圆的周长比直径
eps 浮点相对精度
NaN 空值
NO.2 矩阵运算
----------------------------------
ones 创建一个全为1的矩阵
zeros 创建一个全为0的矩阵
eye 产生单位矩阵
diag(a) 产生对角线元素为a，其余元素为零的矩阵
length 返回矩阵维度的最大值
size 返回矩阵维度
repmat 扩展矩阵
norm 求范数
rand 产生0-1之内的均匀分布的矩阵
randn 产生均值为0，方差为1标准正态分布的矩阵
B.^p 对矩阵B的每一个元素进行操作
randperm 产生范围内的自然数排列
unidrnd 产生范围内随机数，可以指定规模大小
inv 求矩阵的逆
rank 求矩阵的秩
A^p 和 A.^P 对矩阵和矩阵中的元素进行操作
trace 矩阵的迹
matlab 常用快捷键
clc 清屏
clear 从工作空间清除所有变量
cd 改变当前工作子目录
↑： 调用上一次的命令(没有正确执行的命令左边会有一条短横杠)
←： 退后一格
→： 前移一格
who 列出当前工作空间中的变量
whos 列出当前工作空间中的变量及相关信息
Tab 输入几个首字母，按下Tab会自动关联该首字母的所有命令
Ctrl+Tab 在各个窗口切换(没鼠标的时候可能会有点用吧)
（如果自己的程序陷入了死循环怎么办？
如果Matlab一直显示正忙：按下Ctrl+c）
----------------------------------
在代码编辑器中
Ctrl+R 选中多行批量注释
Ctrl+T 取消多行注释
Tab 用来缩进对齐，对多行也有效
F12 在光标所在行设置断点（有鼠标的话直接点左边的短横杠）
F10 单步运行到下一行（我一般直接多设置几个断点，直接点继续）
F5 运行程序
Shift+Enter 可以批量改变量名
Ctrl+C，V 复制粘贴（大家应该都熟悉）
Ctrl+Z 撤销一步
Ctrl+S 保存
Ctrl+A 全选
模拟退火算法
模拟退火算法实例：
问题：请问怎么划分寝室，保证不同寝室间，同寝室同学的共同爱好数量分布尽可能均匀。（方差标准差尽可能小）
共同兴趣爱好： 寝室 N： （A,B）+(A,C)+(A,D)+(B,C)+(B,D)+(C,D)
数据链接：https://wws.lanzous.com/b01htrmwj
python进行数据预处理
import numpy as np
import pandas as pd
df = pd.read_csv(r'load\data.csv') # 读取数据
data = df.iloc[:,1].tolist() # 取出第二列元素，并转化为列表格式
data = [set(i.replace('其他',"").strip('、').split('、'))for i in data] # 将列表中的"其他"替换为空，strip()函数去掉结尾顿号，并且以、为分隔符进行分割并保存到列表中。 最后用set将数据保存到集合中。
result = [len(i&j) for i in data for j in data] #计算每两个集合之间的交集 a&b len(a&b) 交集长度
result = np.array(result).reshape(len(data),-1)
# 转化为numpy中的数组格式，改变形状。 （-1代表交给系统计算）
np.savetxt(r'load\mydata.csv',result,dalimiter = ',')
模拟退火算法
clear ;
clc;
rng('shuffle');% 空置随机数的生成
n = 40; % 学生数量
mydata = csvread('load/mydata.csv');
%% 参数初始化
T0 = 1000; % 初始温度
T = T0; %迭代中
maxgen = 500; %最大迭代个数
LK = 500; %每个温度的迭代速度
alpfa = 0.95; % 温度递减系数
%% 随机生成一个初始解
path0 = randperm(n);
result0 = caculate_ans(path0,mydata); %计算这一组解对应的标准差
%% 定义一些中间变量
min_result = result0;
RESULT = zeros(maxgen,1); % 记录没依次迭代找到的min_result。
%% 模拟退火过程
for iter = 1:maxgen % 外循环，这里采用最大迭代次数
for i = 1:LK % 内循环，在每个温度下开始迭代
path1 = gen_new_path(path0);
result1 = calculate_ans(path1,mydata);
if result1 <result0
path0 = path1;
result0 = result1;
else
p = exp(-(result1 - result0/T));%根据Matroplis准则计算一个概率
if rand(1) < p %生成一个随机数和这个数作比较
path0 = path1;
result0 = result1;
end
end
% 判断是否要更新找到的最佳的解
if result0 < min_result % 如果当前解更好，则对其进行更新
min_result = result0;
best_path = path0;
end
end
RESULT(iter) = min_result; % 保证本轮外循环后找到最优的解
T = alpfa*T; % 温度下降
end
disp('最佳的方案为：');disp(reshape(best_path,4,n/4))
disp('此时的最优值：')；disp(min_result)
%% 画出 迭代产生的图形
figure
plot(1:maxgen,RESULT,'b-');
xlabel('迭代次数');
[~,b] = calculate_ans(best_path,mydata)
toc
function[ans,tmp] = calculate_ans(path,mydata)
n = length(path); %学生个数
nums = n/4;
path = reshape(path,4,nums);
tmp = zeros(1,nums);
for i = 1:nums
tmp(i) = get_total_hobby(path(:,i),mydata);
end
ans = std(tmp); %计算标准差
% ans = (std(tmp) + 1) / max(tmp); 模型改进多目标优化
% ans = (std(tmp) + 1) / mean(tmp);
end
function hobby = get_total_hobby(path,mydata)
% path: 同寝室四名同学的编号
% mydata: 每两个同学之间的共同爱好数
hobby = 0;
for i = 1 : 4
for j = 1:i
if i ~= j
hobby = hobby+mydata(path(i), path(j));
end
end
end
end
function path1 = gen_new_path(path0)
% path0: 原来的路径
n = length(path0);
% 随机选择两种产生新路径的方法
p1 = 0.33; % 使用交换法产生新路径的概率
p2 = 0.33; % 使用移位法产生新路径的概率
r = rand(1); % 随机生成一个[0 1]间均匀分布的随机数
if r< p1 % 使用交换法产生新路径
c1 = randi(n); % 生成1-n中的一个随机整数
c2 = randi(n); % 生成1-n中的一个随机整数
path1 = path0; % 将path0的值赋给path1
path1(c1) = path0(c2); %改变path1第c1个位置的元素为path0第c2个位置的元素
path1(c2) = path0(c1); %改变path1第c2个位置的元素为path0第c1个位置的元素
elseif r < p1+p2 % 使用移位法产生新路径
c1 = randi(n); % 生成1-n中的一个随机整数
c2 = randi(n); % 生成1-n中的一个随机整数
c3 = randi(n); % 生成1-n中的一个随机整数
sort_c = sort([c1 c2 c3]); % 对c1 c2 c3从小到大排序
c1 = sort_c(1); c2 = sort_c(2); c3 = sort_c(3); % c1 < = c2 <= c3
tem1 = path0(1:c1-1);
tem2 = path0(c1:c2);
tem3 = path0(c2+1:c3);
tem4 = path0(c3+1:end);
path1 = [tem1 tem3 tem2 tem4];
else % 使用倒置法产生新路径
c1 = randi(n); % 生成1-n中的一个随机整数
c2 = randi(n); % 生成1-n中的一个随机整数
if c1>c2 % 如果c1比c2大，就交换c1和c2的值
tem = c2;
c2 = c1;
c1 = tem;
end
tem1 = path0(1:c1-1);
tem2 = path0(c1:c2);
tem3 = path0(c2+1:end);
path1 = [tem1 fliplr(tem2) tem3]; %矩阵的左右对称翻转 fliplr，上下对称翻转 flipud
end
end
模拟退火算法原理
https://mp.weixin.qq.com/s/001Klrt7jjf8s5rI3py7Yg
遗传算法
遗传算法实例分析（2）车间调度问题 - 会武术之白猫 - 博客园 (cnblogs.com)
遗传算法解决TSP问题实例
function GA_TCP()
clear all;clc;close all;
%% 1.参数初始化
N=18; %城市的个数
M=200; %种群的个数
C=500; %迭代次数
m=2; %适应值归一化淘汰加速指数
Pc=0.5; %交叉概率
Pmutation=0.2; %变异概率
pos=10*randn(N,2); %生成城市的坐标
D=zeros(N,N); %生成城市之间距离矩阵
%计算两点间的距离
for i=1:N
for j=i+1:N
dis=(pos(i,1)-pos(j,1)).^2+(pos(i,2)-pos(j,2)).^2;
D(i,j)=dis^(0.5);
D(j,i)=D(i,j);
end
end
%% 2.1生成初始群体
popm=zeros(M,N);
for i=1:M
popm(i,:)=randperm(N); %将N序号随机打乱
end
%%2.2随机选择一个种群画出相应的路径
R=popm(1,:);
figure(1);
scatter(pos(:,1),pos(:,2),'ro');
plot_route(pos,R); %画出种群各城市之间的连线
title('18个城市坐标的最初路线图');
str=['初代，','总距离:',num2str(dis)];
text(1000,5000,str);
%% 3.初始化种群及其适应函数
fitness=zeros(M,1);
len=zeros(M,1);%%M为种群的个数
for i=1:M
len(i,1)=myLength(D,popm(i,:));
end
maxlen=max(len);
minlen=min(len);
fitness=fit(len,m,maxlen,minlen);
rr=find(len==minlen);%%在数组len当中找到等于minlen的数组下标记录进rr
R=popm(rr(1,1),:);%%将相应种群当中的编码记录进R
% for i=1:N
% fprintf('%d ',R(i));
% end %%输出对应的编码
% fprintf('\n');
fitness=fitness/sum(fitness);%%适应度函数进行归一化处理
%%这里C为迭代的次数，这里的数组用于记录每次迭代得到的种群中的最小距离
distance_min=zeros(C+1,1); %%用于存储各次迭代的最小的种群的距离
while C>=0
% fprintf('迭代第%d次\n',500-C);
%%选择操作%%
nn=0;
popm_sel=[];
for i=1:size(popm,1)
len_1(i,1)=myLength(D,popm(i,:));
jc=rand*0.8;
for j=1:size(popm,1)
if fitness(j,1)>=jc
nn=nn+1;
popm_sel(nn,:)=popm(j,:);
break;
end
end
end
% %%每次选择都保存最优的种群%%选择函数存在一定的问题，没有达到筛选的效果
% kk=;
% popm_sel=[];
% popm_sel=kk;
[len_m len_index]=min(len_1);
popm_sel=[popm_sel;popm(len_index,:)];%%考虑到随机性，避免最小距离的情况没有被保存，找出最小距离对应的基因添加到popm_sel
%popm_sel(1,:)=popm(len_index,:);
%%交叉操作
nnper=randperm(nn);
for i=1:nn*Pc
A=popm_sel(nnper(i),:);
B=popm_sel(nnper(i+1),:);
[A,B]=cross(A,B);
popm_sel(nnper(i),:)=A;
popm_sel(nnper(i+1),:)=B;
end
%%变异操作
for i=1:nn
pick=rand;
while pick==0
pick=rand;
end
if pick<=Pmutation
popm_sel(i,:)=Mutation(popm_sel(i,:));
end
end
%%求适应度函数
NN=size(popm_sel,1);
len=zeros(NN,1);
for i=1:NN
len(i,1)=myLength(D,popm_sel(i,:));
end
maxlen=max(len);
minlen=min(len);
distance_min(501-C,1)=minlen;
fitness=fit(len,m,maxlen,minlen);
rr=find(len==minlen);%%找到最优基因的下标
%fprintf('minlen=%d\n',minlen);
R=popm_sel(rr(1,1),:);%%找到最优基因
% for i=1:N
% fprintf('%d ',R(i));
% end
% fprintf('\n');
popm=[];
popm=popm_sel;
C=C-1;
%pause(1);
%size(popm)
end
figure(2)
plot_route(pos,R);
title('18个城市坐标的最终优化路线图');
str=['初代，','总距离:',num2str(dis)];
text(1000,5000,str);
figure(3)
plot(distance_min);
end
%% 城市点间连线
function plot_route(a,R)
scatter(a(:,1),a(:,2),'rx');
hold on;
plot([a(R(1),1),a(R(length(R)),1)],[a(R(1),2),a(R(length(R)),2)]);
hold on;
for i=2:length(R)
x0=a(R(i-1),1);
y0=a(R(i-1),2);
x1=a(R(i),1);
y1=a(R(i),2);
xx=[x0,x1];
yy=[y0,y1];
plot(xx,yy);
hold on;
end
end
%% 个体距离计算函数 mylength.m
function len=myLength(D,p) %输入的D为各个点之间的距离矩阵，p为游览顺序的一个一维数组
[N,NN]=size(D);
len=D(p(1,N),p(1,1)); %第一段距离作为初始值
for i=1:(N-1)
len=len+D(p(1,i),p(1,i+1));%%按照顺序把相应的距离相加,相当于对于单个旅行顺序求解距离
end
end
%% 适应度函数fit.m
function fitness=fit(len,m,maxlen,minlen)%%输入的值分别种群当中不同个体对应的距离，适应度加速淘汰指数，最大距离值和最小距离值
fitness=len;
for i =1:length(len)
fitness(i,1)=(1-(len(i,1)-minlen)/(maxlen-minlen+0.0001)).^m;%%为了使得优化速度更快，这里就把目标函数故意求平方
end
end
%% 交叉操作函数 cross.m
function [A,B]=cross(A,B)
L=length(A);
if L<10
W=L;
elseif ((L/10)-floor(L/10))>=rand&&L>10
W=ceil(L/10)+8;
else
W=floor(L/10)+8;
end
%%这里的W为需要交叉的位数，如果小于十则全部交叉，等于10交叉9，大于十则分为不同概率交叉，使得交叉存在较高的概率
p=unidrnd(L-W+1);%%这里是随机产生一个交叉的标志位置
%fprintf('p=%d ',p);
for i=1:W
x=find(A==B(1,p+i-1));
y=find(B==A(1,p+i-1));
[A(1,p+i-1),B(1,p+i-1)]=exchange(A(1,p+i-1),B(1,p+i-1));
[A(1,x),B(1,y)]=exchange(A(1,x),B(1,y));
end
%这里的for循环完成交叉的过程
end
%% 对调函数 exchange.m
function [x,y]=exchange(x,y)
temp=x;
x=y;
y=temp;
end
%% 变异函数 Mutation.m
function a=Mutation(A)
index1=0;
index2=0;
nnper=randperm(size(A,2));
index1=nnper(1);
index2=nnper(2);%%产生两个随即下标，然后选择前两个作为下标进行互换
%fprintf('index1=%d ',index1);
%fprintf('index2=%d ',index2);
temp=0;
temp=A(index1);
A(index1)=A(index2);
A(index2)=temp;
a=A;
end
蚁群算法
%% 蚁群算法学习
clc;
clear all; %清屏清除工作区
close all;
%% 取得要处理的城市的坐标,作出初始图
%position=load('景区经纬度.txt');
%这句可以读入所需要处理的数据，以下采用的是随机产生的坐标信息
t0=clock;
N=30; %设置点数,后面其实有个n和它是一样的，为了理解我都留着
for i = 1 : N
x(i) = rand * 100 ; %生成正态分布的随机坐标点
y(i) = rand * 100 ; %x当作近经度，y当作纬度
scatter(x(i),y(i),'b'); %画出散点图
hold on
text(x(i)+1,y(i)+1,num2str(i)) %用text做好标记，
end
citys=[x;y]';
xlabel('经度'); %横坐标
ylabel('纬度'); %纵坐标
title('TSP蚁群算法优化'); %图片标题
grid on
%% 计算各个城市之间的距离
n = size(citys,1);
D = zeros(n,n);
for i = 1:n
for j = 1:n
if i ~= j
D(i,j) = sqrt(sum((citys(i,:) - citys(j,:)).^2)); %矩阵里元素各自平方后求和
%S(i,j) = sqrt((x(i)-x(j))^2 + (y(i)-y(j))^2); %距离公式
else
D(i,j) = eps; %这里不是零主要是为了之后去的取倒数，对角矩阵值浮点相对精度
end
end
end
%% 种群初始化，参数的设置
m = 75; % 蚂蚁数量
alpha = 1; % 信息素重要程度因子
beta = 5; % 启发函数重要程度因子
vol = 0.2; % 信息素挥发(volatilization)因子
Q = 10; % 常系数
Heu_F = 1./D; % 启发函数(距离的倒数)
Tau = ones(n,n); % 信息素矩阵
Table = zeros(m,n); % 路径记录表（蚂蚁的数量是行数，列数是城市的数量）
%iter = 1; % 迭代次数初值
iter_max = 100; % 最大迭代次数
Route_best = zeros(iter_max,n); % 各代最佳路径（每一次迭代后的最佳路线）
Length_best = zeros(iter_max,1); % 各代最佳路径的长度
Length_ave = zeros(iter_max,1); % 各代路径的平均长度
%% 迭代寻找最佳路径
for iter = 1:iter_max %迭代开始，用for不用再设置加一（原程序是while）
% 随机产生各个蚂蚁的起点城市（原程序）
%start = zeros(m,1);
%for i = 1:m
% temp = randperm(n); %这里产生的是n以内的一个序列
% start(i) = temp(1); %取了这个序列的第一个值，随机值
%end
for i=1:m
Table(i,1) = unidrnd(n); %这个直接产生的就是n以内的随机值
end
% 构建解空间
citys_index = 1:n;
% 逐个蚂蚁选择
for i = 1:m
% 逐个城市路径选择
for j = 2:n
%table是路径记录表，每次j之前的城市都是被选择过的
has_visited = Table(i,1:(j - 1)); % 已访问的城市集合(禁忌表)
allow_index = ~ismember(citys_index,has_visited);
%citys_index中有元素属于禁忌表中元素的时候取1，取反后变成0,产生的是逻辑数组
allow = citys_index(allow_index); % 剩下待访问的城市集合
P = allow;
% 计算城市间转移概率
for k = 1:length(allow)
P(k) = Tau(has_visited(end),allow(k))^alpha * Heu_F(has_visited(end),allow(k))^beta;
%Heu_F是距离的倒数矩阵，所以说越远概率就会越低
end
P = P/sum(P);
% 轮盘赌法选择下一个访问城市
Pc = cumsum(P); %累加函数，把前几个累加到1
target_index = find(Pc >= rand);
target = allow(target_index(1));
%这里和遗传算法不一样,没有采用二分法
Table(i,j) = target;
end
end
%% 首先记录上一次迭代之后的最佳路线放在第一个位置(类似一个学习的效果)
if iter>=2
Table(1,:) = Route_best(iter-1,:);
end
%% 计算各个蚂蚁的路径距离
Length = zeros(m,1);
for i = 1:m
Route = Table(i,:); %取出一条路径
%for循环累加后是第一个点———》最后一个点之间的距离
for j = 1:(n - 1)
Length(i) = Length(i) + D(Route(j),Route(j + 1));
end
%加上最后一个点回到起点的距离
Length(i) = Length(i) + D(Route(n),Route(1));
end
%% 计算最短路径距离及平均距离(这是原程序的代码，感觉冗长。而且我们一般不把判断和选择放在一起)
% if iter == 1
%这部分其实和之前说过的遗传算法best函数的做法类似，先假定第一次产生的就是最佳
%找到最短距离的值的大小和位置
% [min_Length,min_index] = min(Length);
% Length_best(iter) = min_Length; %最佳值先记录下这个最小值
% Length_ave(iter) = mean(Length); %最佳平均值先记录下来这次迭代的平均长度
% Route_best(iter,:) = Table(min_index,:); %记录下这条最优路线
% Limit_iter = 1;
% else
% [min_Length,min_index] = min(Length); %找出最小值和最小值的位置
% Length_best(iter) = min(Length_best(iter - 1),min_Length);
%这个最小值和上个最小值比较，记录最小的那个
% Length_ave(iter) = mean(Length); %求出平均值
%如果下一次迭代值比之前的小，则记录下这条路线，终止迭代次数加一
% if Length_best(iter) == min_Length
% Route_best(iter,:) = Table(min_index,:);
% Limit_iter = Limit_iter+1;
% else
%否则的话路线记录的还是上一条
% Route_best(iter,:) = Route_best((iter-1),:);
% end
% end
[min_Length,min_index] = min(Length); %找到最短距离的值的大小和位置
Length_best(iter,:)= min_Length; %此次迭代的路线最小值记录
Route_best(iter,:) = Table(min_index,:); %此次迭代的路线最小值记录
Length_ave(iter) = mean(Length); %此次迭代的路线平均值记录
%% 更新信息素
Delta_Tau = zeros(n,n);
% 逐个蚂蚁计算
for i = 1:m
% 逐个城市计算
for j = 1:(n - 1)
Delta_Tau(Table(i,j),Table(i,j+1)) = Delta_Tau(Table(i,j),Table(i,j+1)) + Q/Length(i);
%这里可以理解成每只蚂蚁带着等量Q的信息素，均匀撒在路线上
end
%路线最后一个点和起点的信息素
Delta_Tau(Table(i,n),Table(i,1)) = Delta_Tau(Table(i,n),Table(i,1)) + Q/Length(i);
end
Tau = (1-vol)* Tau + Delta_Tau; %信息素挥发一部分再加上增加的
% 迭代次数加1，清空路径记录表
%iter = iter + 1;
Table = zeros(m,n); %清空路线表
end
%% 命令窗口的结果显示
[Shortest_Length,index] = min(Length_best); %找到每次迭代之后记录下来中最小值里面的最小值
Shortest_Route = Route_best(index,:);
Time_Cost=etime(clock,t0); %调用windows系统的时钟进行时间差计算
disp(['环游城市的最短距离:' num2str(Shortest_Length)]);
disp(['最短路径:' num2str([Shortest_Route Shortest_Route(1)])]);
%disp(['收敛迭代次数:' num2str(Limit_iter)]);
disp(['程序执行时间:' num2str(Time_Cost) '秒']);
%% 绘图
plot([citys(Shortest_Route,1);citys(Shortest_Route(1),1)],... %三点省略符为Matlab续行符
[citys(Shortest_Route,2);citys(Shortest_Route(1),2)],'bo-');
%注意最后还要加上起点
text(citys(Shortest_Route(1),1),citys(Shortest_Route(1),2),' 起点');
text(citys(Shortest_Route(end),1),citys(Shortest_Route(end),2),' 终点');
figure(3)
subplot(1,2,1)
plot(Length_best)
xlabel('迭代次数')
ylabel('最短距离收敛轨迹')
subplot(1,2,2)
plot(Length_ave)
xlabel('迭代次数')
ylabel('平均距离收敛轨迹')
鱼群算法
子函数源码：https://mp.weixin.qq.com/s/kHKzKDS6maKuIF7Xrtpt7A
NO.1 觅食行为
function [Xnext,Ynext]=AF_prey(Xi,ii,visual,step,try_number,LBUB,lastY)
%觅食行为
%输入：
%Xi 当前人工鱼的位置
%ii 当前人工鱼的序号
%visual 感知范围
%step 最大移动步长
%try_number 最大尝试次数
%LBUB 各个数的上下限
%lastY 上次的各人工鱼位置的食物浓度
%输出：
%Xnext Xi人工鱼的下一个位置
%Ynext Xi人工鱼的下一个位置的食物浓度
Xnext=[];
Yi=lastY(ii);
for i=1:try_number
Xj=Xi+(2*rand(length(Xi),1)-1)*visual;
%在其视野范围内随机选择一个位置,随机数的范围为（-1，1）
Yj=AF_food(Xj);
if Yi<Yj %求解极大值
Xnext = AF_update(Xi,Xj,step,LBUB);
Xi=Xnext;
break;
end
end
%随机行为
if isempty(Xnext)
Xj=Xi+(2*rand(length(Xi),1)-1)*visual;
Xnext=Xj;
for i=1:length(Xnext)
if Xnext(i)>LBUB(i,2)
Xnext(i)=LBUB(i,2);
end
if Xnext(i)<LBUB(i,1)
Xnext(i)=LBUB(i,1);
end
end
end
Ynext=AF_food(Xnext);
NO.2 聚群行为
----------------------------------
function [Xnext,Ynext]=AF_swarm(X,i,visual,step,deta,try_number,LBUB,lastY)
% 聚群行为
%输入：
%X 所有人工鱼的位置
%i 当前人工鱼的序号
%visual 感知范围
%step 最大移动步长
%deta 拥挤度
%try_number 最大尝试次数
%LBUB 各个数的上下限
%lastY 上次的各人工鱼位置的食物浓度
%输出：
%Xnext Xi人工鱼的下一个位置
%Ynext Xi人工鱼的下一个位置的食物浓度
Xi=X(:,i);
D=dist(Xi',X); %求解欧氏距离，左右矩阵需要列数相等
index=find(D>0 & D<visual);
nf=length(index);
if nf>0 %如果视野内有伙伴
for j=1:size(X,1)
Xc(j,1)=mean(X(j,index)); %求视野内位置点的中
end
Yc=AF_food(Xc);
Yi=lastY(i);
if Yc/nf>deta*Yi
Xnext = AF_update(Xi,Xc,step,LBUB);
Ynext=AF_food(Xnext);
else
[Xnext,Ynext]=AF_prey(Xi,i,visual,step,try_number,LBUB,lastY);
end
else
[Xnext,Ynext]=AF_prey(Xi,i,visual,step,try_number,LBUB,lastY);
end
NO.3 函数值计算
----------------------------------
function Y = AF_food(X)
%计算人工鱼当前位置的食物浓度,即为要处理的函数表达式
num = size(X,2);
Y = zeros(1,num);
for i=1:num
Y(1,i) = (sin(X(1,i))/X(1,i))*(sin(X(2,i))/X(2,i));
end
NO.4 位置更新
----------------------------------
function Xnext = AF_update(Xi,Xcho,step,LBUB)
%用来更新位置
Xnext=Xi+rand*step*(Xcho-Xi)/norm(Xcho-Xi); %norm求解范数，默认二范数
for i=1:length(Xnext)
if Xnext(i)>LBUB(i,2)
Xnext(i)=LBUB(i,2);
end %对边界的处理，越界则拉回
if Xnext(i)<LBUB(i,1)
Xnext(i)=LBUB(i,1);
end
end
NO.5 初始化函数
----------------------------------
function X=AF_init(Nfish,lb_ub)
%输入：
% Nfish 鱼群大小
% lb_ub 鱼的活动范围
%输出：
% X 产生的初始人工鱼群
%%lb_ub每行中前两个数是范围的上下限，第3个数自变量的维数
lb=lb_ub(1);
ub=lb_ub(2);
nr=lb_ub(3);
X=zeros(nr,Nfish);
for j=1:nr
X(j,:)=lb+(ub-lb)*rand(1,Nfish);
%X(j,:) = unifrnd(lb,ub,[1,Nfish]); 其实直接产生范围内的随机数也可以
end
NO.6 追尾函数
----------------------------------
function [Xnext,Ynext]=AF_follow(X,i,visual,step,deta,try_number,LBUB,lastY)
% 追尾行为
%输入：
%X 所有人工鱼的位置
%i 当前人工鱼的序号
%visual 感知范围
%step 最大移动步长
%deta 拥挤度
%try_number 最大尝试次数
%LBUB 各个数的上下限
%lastY 上次的各人工鱼位置的食物浓度
%输出：
%Xnext Xi人工鱼的下一个位置
%Ynext Xi人工鱼的下一个位置的食物浓度
Xi=X(:,i);
D=dist(Xi',X);
index=find(D>0 & D<visual);
nf=length(index);
if nf>0
XX=X(:,index);
YY=lastY(index);
[Ymax,Max_index]=max(YY);
Xmax=XX(:,Max_index);
Yi=lastY(i);
if Ymax/nf>deta*Yi;
Xnext = AF_update(Xi,Xmax,step,LBUB);
Ynext=AF_food(Xnext);
else
[Xnext,Ynext]=AF_prey(X(:,i),i,visual,step,try_number,LBUB,lastY);
end
else
[Xnext,Ynext]=AF_prey(X(:,i),i,visual,step,try_number,LBUB,lastY);
end
用来对比的一个小程序
----------------------------------
clc
clear all;
[X,Y] = meshgrid(-10:0.1:10);
Z=(sin(X)./X).*(sin(Y)./Y);
mesh(X,Y,Z)
hold on
[max_value,index] = max(Z(:))
[I_row, I_col] = ind2sub(size(Z),index);
%ind2sub将线性化索引转换成下标
disp(['最大值点X：',num2str(X(I_row,I_col))]); %因为函数是对称的
disp(['最大值点Y：',num2str(Y(I_row,I_col))]);
plot3(X(I_row,I_col),Y(I_row,I_col),max_value,'r*','MarkerSize',10)
%plot3(-0.0068486,-0.0068486,0.99998,'bo','MarkerSize',10)
%可以把鱼群算法的最大值点数值代入对比一下
主函数源码： https://mp.weixin.qq.com/s/MDqdIuLqCGBJAQ-1bcnJwg
主函数
----------------------------------
clc
clear all
close all
tic
%% 参数设置
fishnum=100; %生成100只人工鱼
MAXGEN=100; %最多迭代次数
try_number=100; %最多试探次数
visual=1; %感知距离
delta=0.618; %拥挤度因子
step=0.1; %步长
%% 初始化鱼群
lb_ub=[-10,10,2];
X=AF_init(fishnum,lb_ub);
LBUB=repmat(lb_ub(1:2),lb_ub(1,3),1); %用作边界判断
gen=1;
times=1;
BestY=zeros(1,MAXGEN); %每步中最优的函数值
BestX=zeros(2,MAXGEN); %每步中最优的自变量
besty=-inf; %最优函数值
Y=AF_food(X);
%% 初始状态图
[X_plot,Y_plot] = meshgrid(-10:0.1:10);
Z=(sin(X_plot)./X_plot).*(sin(Y_plot)./Y_plot);
mesh(X_plot,Y_plot,Z);
xlabel('第一维度x的取值范围');
ylabel('第二维度y的取值范围');
zlabel('函数值');
hold on;
%% 算法开始
while gen<=MAXGEN
for i=1:fishnum %每个坐标点都更新一次
[Xi1,Yi1]=AF_swarm(X,i,visual,step,delta,try_number,LBUB,Y); %聚群行为
[Xi2,Yi2]=AF_follow(X,i,visual,step,delta,try_number,LBUB,Y); %追尾行为
%这两种行为的Xnext更新方式类似，保留一个就好
if Yi1>Yi2 %选择聚群和追尾行为的较大值保留
X(:,i)=Xi1;
Y(1,i)=Yi1;
else
X(:,i)=Xi2;
Y(1,i)=Yi2;
end
end
if times==2 %两次刷新一次
cla;
mesh(X_plot,Y_plot,Z);
plot3(X(1,:),X(2,:),Y,'o','color',[gen/MAXGEN,0,0])
title('鱼群算法动态搜寻');
pause(0.2);
times=0;
end
%颜色根据迭代次数而颜色不同
[Ymax,index]=max(Y);
if Ymax>besty
besty=Ymax; %更新一下最大值
bestx=X(:,index);
BestY(gen)=Ymax;
BestX(:,gen)=X(:,index);
else
BestY(gen)=BestY(gen-1); %如果没有上次值大则保持
BestX(:,gen)=BestX(:,gen-1);
end
gen=gen+1;
times=times+1;
end
figure
mesh(X_plot,Y_plot,Z);
hold on;
plot3(bestx(1),bestx(2),besty,'r*','MarkerSize',10)
title('鱼群算法最大值点')
figure
plot(BestY)
xlabel('迭代次数')
ylabel('优化值')
title('鱼群算法迭代过程')
disp(['最大值位置：',num2str(bestx')])
disp(['最大函数值：',num2str(besty,'%1.5f')])
toc;
粒子群算法
NO.1 二维PSO算法代码
NO.1 二维PSO算法代码
clc;clear;close all; %清屏，删除已有的变量和窗口,tic和toc用来给程序计时
%% 函数的定义；
%f= @(x)x .* sin(x) .* cos(2 * x) - 2 * x.* sin(3 * x); % 函数表达式
%figure(1);ezplot(f,[0,20]); %ezplot只用来画一次函数
f= @(a,b)(a .* sin(a) .* cos(2 * a) - 2 * a .* sin(3 * a)).*(b .* sin(b) .* cos(2 * b) - 2 * b .* sin(3 * b)); % 函数表达式
figure(1);
[x0_1, x0_2]=meshgrid(0:.2:20); %生成网格数据，即为一些采样点
y0=f(x0_1,x0_2); %将这些值带入函数计算
C=gradient(y0);
mesh(x0_1, x0_2,y0,C);
% colorbar;%可以看到色板，也可以给指定区域指定颜色
%set(h,'EdgeColor','b','FaceColor','w','MarkerEdgecolor','r','MarkerFacecolor','y')
xlabel('第一维度x的取值范围');
ylabel('第二维度y的取值范围');
zlabel('函数值');
hold on;
%% 开始种群等基本定义
N = 500; % 初始种群个数
d = 2; % 空间维数（参看上述的函数表达式）
ger = 300; % 最大迭代次数
plimit = [0,20;0,20]; % 设置位置参数限制(矩阵的形式可以多维)，现在2X2矩阵
vlimit = [-1.5, 1.5;-1.5, 1.5]; % 设置速度限制
w = 0.8; % 惯性权重,个体历史成绩对现在的影响0.5~1之间
%还有自适应调整权重、随机权重等等
%（不同的权重设置很影响性能，按需要选取）
c1 = 0.5; % 自我学习因子
c2 = 0.5; % 群体学习因子
tic; %计时开始
for i = 1:d
x(:,i) = plimit(i, 1) + (plimit(i, 2) - plimit(i, 1)) * rand(N, 1);%初始种群的位置
end %rand(N,1)产生N行一列范围在1之内的随机数
%第一列,第二列：x=0+（20-0）*（1之内的随机数）
v = rand(N, d); % 初始种群的速度,500行2列分别在两个维度上
xm = x; % 每个个体的历史最佳位置
ym = zeros(1, d); % 种群的历史最佳位置，两个维度，设置为0
fxm = zeros(N, 1); % 每个个体的历史最佳适应度，设置为0
fym = -inf; % 种群历史最佳适应度,求最大值先设置成负无穷
plot3(xm(:,1),xm(:,2),f(xm(:,1),xm(:,2)), 'mo'); %r红,g绿,b蓝,c蓝绿,m紫红,y黄,k黑,w白
title('种群初始分布状态图'); %plot3在三维区域画出空间上的点，把格式设置成‘o’用来画每个位置的散点图
hold on;
%figure(2);
%mesh(x0_1, x0_2, y0);
%hold on;
%plot3(xm(:,1),xm(:,2),f(xm(:,1),xm(:,2)), 'ro');
%hold on;
iter=1; %初始的迭代次数因为用while设置为一
times = 1;
record = zeros(ger, 1); %记录器
%% 迭代更新开始
while iter <= ger
fx = f(x(:,1),x(:,2)); % 代入x中的二维数据，算出个体当前适应度,为500行1列的数据
for i = 1:N %对每一个个体做判断
if fxm(i) < fx(i) %如果每个个体的历史最佳适应度小于个体当前适应度
fxm(i) = fx(i); % 更新个体历史最佳适应度,第一轮就是把小于零的清除
xm(i,:) = x(i,:); % 更新个体历史最佳位置
end
end
if fym < max(fxm) %种群历史最佳适应度小于个体里面最佳适应度的最大值
[fym, nmax] = max(fxm); % 更新群体历史最佳适应度,取出最大适应度的值和所在行数即位置
ym = xm(nmax, :); % 更新群体历史最佳位置
end
v = v * w + c1 * rand *(xm - x) + c2 * rand *(repmat(ym, N, 1) - x); % 速度更新公式,repmat函数把ym矩阵扩充成N行1列
%%边界速度处理
for i=1:d
for j=1:N
if v(j,i)>vlimit(i,2); %如果速度大于边界速度，则把速度拉回边界
v(j,i)=vlimit(i,2);
end
if v(j,i) < vlimit(i,1) %如果速度小于边界速度，则把速度拉回边界
v(j,i)=vlimit(i,1);
end
end
end
x = x + v; % 位置更新
for i=1:d
for j=1:N
if x(j,i)>plimit(i,2)
x(j,i)=plimit(i,2);
end
if x(j,i) < plimit(i,1)
x(j,i)=plimit(i,1);
end
end
end
record(iter) = fym; %记录最大值
if times >= 10
cla; %清除轴线图形
mesh(x0_1, x0_2, y0);
plot3(x(:,1),x(:,2),f(x(:,1),x(:,2)), 'ro');title('状态位置变化');
pause(0.5);
times=0;
end
iter = iter+1;
times=times+1;
end
%% 作图
figure(3);
plot(record); %画出最大值的变化过程
title('收敛过程');
figure(4);
mesh(x0_1, x0_2, y0);
hold on;
plot3(x(:,1),x(:,2),f(x(:,1),x(:,2)), 'ro');title('最终状态图');
disp(['最大值为：',num2str(fym)]);
disp(['最大值点位置：',num2str(ym)]);
toc; %计时结束
%% 关于权重问题参考CSDN博主
%https://blog.csdn.net/qq_41053605/article/details/89156125
----------------------------------
NO.2 三维PSO
NO.2 三维PSO
%改动后的3维空间的粒子群算法2019/08/11
clc;clear;close all;
%% 初始化种群
f= @(a,b,c) a.* sin(a) + b.^2 + 3 * c; % 函数表达式
%figure(1);
%[x0_1, x0_2]=meshgrid(0:.2:20);
%y0=f(x0_1,x0_2);
%mesh(x0_1, x0_2, y0);
%hold on;
N = 500; % 初始种群个数
d = 3; % 空间维数
ger = 300; % 最大迭代次数
w=0.9;
limit = [-2, 2;-2,2; -2,2]; % 设置位置参数限制(矩阵的形式可以多维)
vlimit = [-1.5, 1.5;-1.5, 1.5; -1.5, 1.5]; % 设置速度限制
% 惯性权重
c_2 = 0.5; % 自我学习因子
c_3 = 0.5; % 群体学习因子
for i = 1:d
x(:,i) = limit(i, 1) + (limit(i, 2) - limit(i, 1)) * rand(N, 1);%初始种群的位置
end
v = rand(N, d); % 初始种群的速度
xm = x; % 每个个体的历史最佳位置
ym = zeros(1, d); % 种群的历史最佳位置
fxm = zeros(N, 1); % 每个个体的历史最佳适应度
fym = -inf; % 种群历史最佳适应度
%plot3(xm(:,1),xm(:,2),f(xm(:,1),xm(:,2)), 'ro');title('初始状态图');
%hold on;
%figure(2);
%mesh(x0_1, x0_2, y0);
%hold on;
%plot3(xm(:,1),xm(:,2),f(xm(:,1),xm(:,2)), 'ro');
%hold on;
%% 粒子群工作
tic;
iter = 1;
times = 1;
record = zeros(ger, 1); % 记录器
while iter <= ger
fx = f(x(:,1),x(:,2),x(:,3)) ; % 个体当前适应度
for i = 1:N
if fxm(i) < fx(i)
fxm(i) = fx(i); % 更新个体历史最佳适应度
xm(i,:) = x(i,:); % 更新个体历史最佳位置
end
end
if fym < max(fxm)
[fym, nmax] = max(fxm); % 更新群体历史最佳适应度
ym = xm(nmax, :); % 更新群体历史最佳位置
end
v = v * w + c_2 * rand *(xm - x) + c_3 * rand *(repmat(ym, N, 1) - x);% 速度更新
% 边界速度处理
for i=1:d
for j=1:N
if v(j,i)>vlimit(i,2)
v(j,i)=vlimit(i,2);
end
if v(j,i) < vlimit(i,1)
v(j,i)=vlimit(i,1);
end
end
end
x = x + v;% 位置更新
% 边界位置处理
for i=1:d
for j=1:N
if x(j,i)>limit(i,2)
x(j,i)=limit(i,2);
end
if x(j,i) < limit(i,1)
x(j,i)=limit(i,1);
end
end
end
record(iter) = fym;%最大值记录
% if times >= 10
% cla;
% mesh(x0_1, x0_2, y0);
% plot3(x(:,1),x(:,2),f(x(:,1),x(:,2)), 'ro');title('状态位置变化');
% pause(0.5);
%times=0;
%end
iter = iter+1;
times=times+1;
end
%figure(3);plot(record);title('收敛过程')
%figure(4);
%mesh(x0_1, x0_2, y0);
%hold on;
%plot3(x(:,1),x(:,2),f(x(:,1),x(:,2)), 'ro');title('最终状态图');
disp(['最大值：',num2str(fym)]);
disp(['变量取值：',num2str(ym)]);
toc;
遗传算法国赛使用案例：
%主函数
function main()
clear;clc;
popsize=100;%种群大小
chromlength=15;%二进制码长 1—6 质量比 7—15 链长
pc=0.6;%交叉概率
pm=0.005;%变异概率32
rou=3.2;%链密度
v=32.6;%风速
V=0.55;%水速
beta=0;%风速水速夹角
%H=20;%水深
pop=initpop(popsize,chromlength,rou,v,V,beta);%初始种群
%-----------------------------------------------------------------------------
for i=1:200%最大遗传代数 200
[objvalue]=cal_objvalue(pop,v,V,beta,rou);%计算适应度值（ 函数值）
fitvalue=objvalue;
%size(pop)
[newpop]=selection(pop,fitvalue);%选择优秀个体产生新种群
%size(newpop)
[newpop]=crossover(newpop,pc,rou,v,V,beta);%对新种群进行概率为 PC 的交叉操作
[newpop]=mutation(newpop,pm,rou,v,V,beta);%对新种群进行概率为 pm 的变异操作
pop=newpop;%更新种群
[bestindividual,~]=best(pop,fitvalue);%寻找最优解
bestr=binary2decimal(bestindividual(1,1:6));%binary2decimal 将质量比转换为十进制
bestL=binary2decimal(bestindividual(1,7:15));%binary2decimal 将链长转换为十进制
end
[phi0,~,~,~,~,h,theta,x]=moor(bestr,20,v,V,beta,rou,bestL);
fprintf(' 锚链密度是-->%5.2f\n',rou);
fprintf(' 铁球 r-->%5.2f\n',bestr);
fprintf(' 锚链长 L-->%5.2f\n',bestL);
fprintf(' 吃水 h-->%5.2f\n',h);
fprintf(' 浮动范围 x-->%5.2f\n',x);
fprintf(' 铁通倾角 phi0-->%5.2f\n',rad2deg(phi0));
fprintf(' 海床夹角 theta-->%5.2f\n',rad2deg(theta));
fprintf(' 水速 V-->%5.2f\n',V);
fprintf(' 风速 v-->%5.2f\n',v);
fprintf(' 水速风速夹角 beta-->%5.2f\n',beta);
end
%选择优秀个体产生新种群用轮盘赌法选择
function [newpop]=selection(pop,fitvalue)
[px,~]=size(pop);
totalfit=sum(fitvalue);
p_fitvalue=fitvalue/totalfit;%用每个个体的适应度除以总适应度使个体的概率在 0-1
p_fitvalue=cumsum(p_fitvalue);%一列累计求和
ms=sort(rand(px,1));%随机生成 100*1 的随机数后排大小
fitin=1;
newin=1;
p_fitvalue(100)=1;33
while newin<=px %比较 1-100
if ms(newin)<=p_fitvalue(fitin)%个体概率和大于随机数
newpop(newin,:)=pop(fitin,:);%将满足条件的优秀个体加入到新种群中
newin=newin+1;
else fitin=fitin+1;
end
end
End
%寻找最优个体和最优适应度 y 值
function [bestindividual,bestfit]=best(pop,fitvalue)
[px,~]=size(pop);
bestindividual=pop(1,:);
bestfit=fitvalue(1);
for i=2:px
if fitvalue(i)>bestfit%更新最优个体
bestfit=fitvalue(i);
bestindividual=pop(i,:);
end
end
End
%二进制转化为十进制
function pop2=binary2decimal(pop)
[~,py]=size(pop);
for i=1:py
pop1(:,i)=2.^(py-i).*pop(:,i);
end
temp=sum(pop1,2);
if py==6
pop2=0.5+temp*4.5/2^6;%质量比 r 0.5-5
elseif py==9
pop2=10+temp*30/2^9;%链长度 L 10-40
end
end
%计算适应度
function [objvalue]=cal_objvalue(pop,v,V,beta,rou)
r=binary2decimal(pop(:,1:6));%质量比解码
L=binary2decimal(pop(:,7:15));%链长度解码
%size(r)
objvalue=ones(100,1);
for i=1:10034
[phi0,~,~,~,~,h,~,x]=moor(r(i),20,v,V,beta,rou,L(i));%带入函数求值
objvalue(i,:)=20-0.1701*x-0.8299*h-0*phi0;%钢桶倾角 phi0,游动区域 x,吃水深度 h
end
end
%对挑选个体进行交叉
function [newpop]=crossover(pop,pc,rou,v,V,beta)
[px,py]=size(pop);
newpop=ones(size(pop));
for i=1:2:px-1
if(rand<pc)%比较交叉概率
phi0_1=1;
theta_1=1;
phi0_2=1;
theta_2=1;
while((rad2deg(phi0_1)>5)||(rad2deg(phi0_2)>5)||(rad2deg(theta_1)>16)||(rad2deg(theta_2)>16)
)
choice=round(rand);%choice 为 0 在前部交叉， choice 为 1 在后部交叉
if(choice==1)
cpoint=round(rand*6);%随机产生交叉位置
newpop(i,:)=[pop(i,1:cpoint),pop(i+1,cpoint+1:py)];%基因交换
newpop(i+1,:)=[pop(i+1,1:cpoint),pop(i,cpoint+1:py)];%基因交换
else
cpoint=6+round(rand*9);%随机产生交叉位置
newpop(i,:)=[pop(i,1:cpoint),pop(i+1,cpoint+1:py)];%基因交换
newpop(i+1,:)=[pop(i+1,1:cpoint),pop(i,cpoint+1:py)];%基因交换
end
r_1=binary2decimal(newpop(i,1:6));%质量比解码
L_1=binary2decimal(newpop(i,7:15));%链长度解码
[phi0_1,~,~,~,~,~,theta_1,~]=moor(r_1,20,v,V,beta,rou,L_1);%带入函数求值
r_2=binary2decimal(newpop(i+1,1:6));%质量比解码
L_2=binary2decimal(newpop(i+1,7:15));%链长度解码
[phi0_2,~,~,~,~,~,theta_2,~]=moor(r_2,20,v,V,beta,rou,L_2);%带入函数求值
end
else %不变异保持原来基因
newpop(i,:)=pop(i,:);
newpop(i+1,:)=pop(i+1,:);
end
end
end35
%随机生成 100 个 10 位二进制代码组初始化种群
function pop=initpop(popsize,chromlength,rou,v,V,beta)%popsize=100,chromlength=16
pop=ones(100,chromlength);
for i=1:popsize
a=round(rand(1,chromlength));%rand 随机生成 0-1 的随机数 round 四舍五入生成 01
r=binary2decimal(a(1,1:6));%质量比解码
L=binary2decimal(a(1,7:15));%链长度解码
[phi0,~,~,~,~,~,theta,~]=moor(r,20,v,V,beta,rou,L);%带入函数求值
if((rad2deg(phi0)<=5)&&(rad2deg(theta)<=16))
[phi0,~,~,~,~,~,theta,~]=moor(r,16,v,V,beta,rou,L);%带入函数求值
if((rad2deg(phi0)<=5)&&(rad2deg(theta)<=16))
pop(i,:)=a(1,:);%满足角度条件则加入种群
else
i=i-1;
end
end
end
%size(pop)
end
%挑选个体进行变异
function [newpop]=mutation(pop,pm,rou,v,V,beta)
[px,py]=size(pop);
newpop=ones(size(pop));
for i=1:px
phi0=1;
theta=1;
if(rand<pm)%比较变异概率
while((rad2deg(phi0)>5)||(rad2deg(theta)>16))%检验变异后是否满足约束条件若不满足重新进
行变异
mpoint=round(rand*py);%找到变异基因位置
if mpoint<=0
mpoint=1;
end
newpop(i,:)=pop(i,:);
if newpop(i,mpoint)==0
newpop(i,mpoint)=1;
else newpop(i,mpoint)=0;
end
r=binary2decimal(newpop(i,1:6));%质量比解码
L=binary2decimal(newpop(i,7:15));%链长度解码36
[phi0,~,~,~,~,~,theta,~]=moor(r,20,v,V,beta,rou,L);%带入函数求值
end
else %没有变异
newpop(i,:)=pop(i,:);
end
end
end
1:6));%质量比解码
L_2=binary2decimal(newpop(i+1,7:15));%链长度解码
[phi0_2,~,~,~,~,~,theta_2,~]=moor(r_2,20,v,V,beta,rou,L_2);%带入函数求值
end
else %不变异保持原来基因
newpop(i,:)=pop(i,:);
newpop(i+1,:)=pop(i+1,:);
end
end
end35
%随机生成 100 个 10 位二进制代码组初始化种群
function pop=initpop(popsize,chromlength,rou,v,V,beta)%popsize=100,chromlength=16
pop=ones(100,chromlength);
for i=1:popsize
a=round(rand(1,chromlength));%rand 随机生成 0-1 的随机数 round 四舍五入生成 01
r=binary2decimal(a(1,1:6));%质量比解码
L=binary2decimal(a(1,7:15));%链长度解码
[phi0,~,~,~,~,~,theta,~]=moor(r,20,v,V,beta,rou,L);%带入函数求值
if((rad2deg(phi0)<=5)&&(rad2deg(theta)<=16))
[phi0,~,~,~,~,~,theta,~]=moor(r,16,v,V,beta,rou,L);%带入函数求值
if((rad2deg(phi0)<=5)&&(rad2deg(theta)<=16))
pop(i,:)=a(1,:);%满足角度条件则加入种群
else
i=i-1;
end
end
end
%size(pop)
end
%挑选个体进行变异
function [newpop]=mutation(pop,pm,rou,v,V,beta)
[px,py]=size(pop);
newpop=ones(size(pop));
for i=1:px
phi0=1;
theta=1;
if(rand<pm)%比较变异概率
while((rad2deg(phi0)>5)||(rad2deg(theta)>16))%检验变异后是否满足约束条件若不满足重新进
行变异
mpoint=round(rand*py);%找到变异基因位置
if mpoint<=0
mpoint=1;
end
newpop(i,:)=pop(i,:);
if newpop(i,mpoint)==0
newpop(i,mpoint)=1;
else newpop(i,mpoint)=0;
end
r=binary2decimal(newpop(i,1:6));%质量比解码
L=binary2decimal(newpop(i,7:15));%链长度解码36
[phi0,~,~,~,~,~,theta,~]=moor(r,20,v,V,beta,rou,L);%带入函数求值
end
else %没有变异
newpop(i,:)=pop(i,:);
end
end
end
4130
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
