选择操作：从旧个体中以一定概率选择优良个体组成新的种群，以繁殖得到下一代。
轮盘赌法：即基于适应度比例的选择策略，个体i被选中的概率为
交叉操作：从种群中随机选择两个个体，通过两个染色体的交换组合，把父串的优秀特征遗传给子串，从而产生新的优秀个体。
采用实数交叉，第k个染色体ak和第l个染色体al在j位的交叉操作方法为,b为[0, 1]随机数：
变异操作：从种群中随机选择一个个体，选择个体中的一点进行变异以产生更优秀的个体。
第i个个体的第j个基因aij进行变异操作的方法为,r为[0,1]随机数,gen为当前迭代次数，genmax为最大迭代次数：
遗传算法代码很多，并且运算速度慢，相比于下面的粒子群算法更复杂，所以如果遇到优化问题，粒子群算法更方便。
核心在于更新粒子速度和位置
算法流程：
由n只鸟组成的种群
B
=
(
B
1
,
B
2
,
B
3
,
…
,
B
n
)
B = (B_1, B_2, B_3,…, B_n)
B=(B1,B2,B3,…,Bn)，其中第
i
i
i只鸟
B
i
Bi
Bi的位置可表示为
X
i
=
[
x
i
1
,
x
i
2
,
…
,
x
i
d
]
X_i = [x_{i1}, x_{i2}, …, x_{id}]
Xi=[xi1,xi2,…,xid];其速度为
V
i
=
[
v
i
1
,
v
i
2
,
…
,
V
i
d
]
Vi = [v_{i1}, v_{i2}, …, V_{id}]
Vi=[vi1,vi2,…,Vid]，则第
t
+
1
t+1
t+1次迭代中鸟类的位置与速度可更新为：
新速度= 惯性项（w V）+个体自信度（i best）+全局最优（best）
W 一般是1或者0.98，也可以随着鸟的飞行不断变小
c1 = c2 = 2或者0.5
fun函数：
function y = fun(x)
y = sum(x .^ 2);
if (x(1) + x(2) > 2) && (x(1) - x(2) > 0.5)
y = y;
else
y = y + 1000;
end
粒子群算法：
% 清空
clc
clear
%% 1、初始化参数
c1 = 1.49445;
c2 = 1.49445;
w = 1;
maxgen = 200; % 迭代次数
sizepop = 50; % 粒子个数
nVar = 2;
% 速度上下限
Vmax = 5;
Vmin = -5;
% 位置上下限
popmax = 100;
popmin = -100;
%% 2、初始化粒子位置
for i = 1 : sizepop
pop(i, :) = (popmax - popmin) * rands(1, nVar) + popmin;
V(i, :) = (Vmax - Vmin) * rands(1, nVar) + Vmin;
fitness(i) = fun(pop(i, :)); % 3、计算粒子目标函数值
end
%% 4、确定全局最优粒子
[bestfitness, bestindex] = min(fitness);
gbest = pop(bestindex, :);
pbest = pop;
fitnesspbest = fitness;
fitnessgbest = bestfitness;
%% 5、更新粒子速度和位置
for i = 1 : maxgen
for j = 1 : sizepop
w = 0.98 .^ i * w;
% 速度更新
V(j, :) = w .* V(j, :) + c1 * rand * (pbest(j, :) - pop(j, :)) + c2 * rand * (gbest - pop(j, :));
% 限定速度大小
V(j, find(V(j, :) > Vmax)) = Vmax;
V(j, find(V(j, :) < Vmin)) = Vmin;
% 更新位置
pop(j, :) = pop(j, :) + V(j, :);
% 限定位置
pop(j, find(pop(j, :) > popmax)) = popmax;
pop(j, find(pop(j, :) < popmin)) = popmin;
% 6、计算粒子目标函数值
fitness(j) = fun(pop(j, :));
end
% 更新全局最优粒子
for j = 1 : sizepop
if fitness(j) < fitnesspbest(j)
pbest(j, :) = pop(j, :);
fitnesspbest(j) = fitness(j);
end
if fitness(j) < fitnessgbest
gbest = pop(j, :);
fitnessgbest = fitness(j);
end
end
yy(i) = fitnessgbest;
end
%% 画图
plot(yy)
title('最优个体适应度', 'fontsize', 12);
xlabel('进化代数', 'fontsize', 12);
ylabel('适应度', 'fontsize', 12);
disp('函数值 变量');
disp([fitnessgbest, gbest]);
函数值 变量
2.1251 1.2501 0.7499
惩罚函数：用于解决约束条件下的最优化问题。通过惩罚函数可以将有约束的目标函数转化为无约束的目标函数。
f
(
X
)
′
=
f
(
X
)
+
C
f(X)^{'} = f(X)+C
f(X)′=f(X)+C
在实际问题中大都具有多个目标且需要同时满足，即在同一问题模型中同时存在几个非线性目标，而这些目标函数需要同时进行优化处理，并且这些目标又往往是互相冲突的，称这类问题称为多目标规划问题。
支配：
在多目标优化问题中，如果个体p至少有一个目标比个体q好，而且个体p的所有目标都不比个体q差，那么称个体p支配个体q。
A和B支配C，但是A和B不存在支配关系
序值：如果p支配q，那么p的序值比q低。如果p和q互不支配，那么p和q有相同的序值。
拥挤距离：用来计算某前端中的某个体与该前端中其他个体之间的距离，用以表征个体间的拥挤程度。
example.m
clear
clc
MultiObj.fun = @(x) [-10.*(exp(-0.2.*sqrt(x(:,1).^2+x(:,2).^2)) + exp(-0.2.*sqrt(x(:,2).^2+x(:,3).^2))), ...
sum(abs(x).^0.8 + 5.*sin(x.^3),2)];
MultiObj.nVar = 3;
MultiObj.var_min = -5 .* ones(1, MultiObj.nVar);
MultiObj.var_max = 5 .* ones(1, MultiObj.nVar);
params.Np = 200;
params.Nr = 150; % 想要多少个点
params.maxgen = 100;
params.W = 0.4;
params.C1 = 2;
params.C2 = 2;
params.ngrid = 20; % 网格个数。算拥挤距离
params.maxvel = 5;
params.u_mut = 0.5; % 变异率（防止粒子群算法陷入局部最优）
REP = MOPSO(params, MultiObj);
MOPSO.m
function REP = MOPSO(params, MultiObj)
Np = params.Np;
Nr = params.Nr;
maxgen = params.maxgen;
W = params.W;
C1 = params.C1;
C2 = params.C2;
ngrid = params.ngrid;
maxvel = params.maxvel;
u_mut = params.u_mut;
fun = MultiObj.fun; % 优化函数
nVar = MultiObj.nVar; % 变量个数
var_min = MultiObj.var_min(:);
var_max = MultiObj.var_max(:);
POS = repmat((var_max-var_min)', Np, 1) .* rand(Np, nVar) + repmat(var_min', Np, 1);
VEL = zeros(Np, nVar); % 初始速度为0
POS_fit = fun(POS); % 计算适应度
PBEST = POS;
PBEST_fit = POS_fit;
DOMINATED = checkDomination(POS_fit);
REP.pos = POS(~DOMINATED, :);
REP.pos_fit = POS_fit(~DOMINATED, :);
REP = updateGrid(REP, ngrid);
maxvel = (var_max - var_min) .* maxvel ./ 100;
gen = 1;
display(['Generation #0 - Repository size: ' num2str(size(REP.pos,1))]);
stopCondition = false;
while ~stopCondition
h = selectLeader(REP);
VEL = W .* VEL + C1 * rand(Np, nVar) .* (PBEST - POS) ...
+ C2 * rand(Np, nVar) .* (repmat(REP.pos(h, :), Np, 1) - POS);
POS = POS + VEL;
POS = mutation(POS, gen, maxgen, Np, var_max, var_min, nVar, u_mut);
[POS, VEL] = checkBoundaries(POS, VEL, maxvel, var_max, var_min);
POS_fit = fun(POS);
REP = updateRepository(REP, POS, POS_fit, ngrid); % 更新记忆库
if(size(REP.pos, 1) > Nr) % 基因库大于所需个体
REP = deleteFromRepository(REP, size(REP.pos, 1) - Nr, ngrid);
end
pos_best = dominates(POS_fit, PBEST_fit);
best_pos = ~dominates(PBEST_fit, POS_fit);
best_pos(rand(Np, 1) >= 0.5) = 0;
if(sum(pos_best) > 1)
PBEST_fit(pos_best, :) = POS_fit(pos_best, :);
PBEST(pos_best, :) = POS(pos_best, :);
end
if(sum(best_pos) > 1)
PBEST_fit(best_pos, :) = POS_fit(best_pos, :);
PBEST(best_pos, :) = POS(best_pos, :);
end
display(['Generation #' num2str(gen) ' - Repository size: ' num2str(size(REP.pos,1))]);
gen = gen + 1;
if(gen > maxgen)
stopCondition = true;
end
end
h_rep = plot(REP.pos_fit(:,1),REP.pos_fit(:,2),'ok');
end
function REP = updateRepository(REP, POS, POS_fit, ngrid)
% 更新基因库
DOMINATED = checkDomination(POS_fit);
REP.pos = [REP.pos; POS(~DOMINATED, :)];
REP.pos_fit = [REP.pos_fit; POS_fit(~DOMINATED, :)];
DOMINATED = checkDomination(REP.pos_fit);
REP.pos_fit = REP.pos_fit(~DOMINATED, :);
REP.pos = REP.pos(~DOMINATED, :);
REP = updateGrid(REP, ngrid); % 更新网格
end
function [POS, VEL] = checkBoundaries(POS, VEL, maxvel, var_max, var_min)
% 边界控制
Np = size(POS, 1);
MAXLIM = repmat(var_max(:)', Np, 1);
MINLIM = repmat(var_min(:)', Np, 1);
MAXVEL = repmat(maxvel(:)', Np, 1);
MINVEL = repmat(-maxvel(:)', Np, 1);
VEL(VEL > MAXVEL) = MAXVEL(VEL > MAXVEL);
VEL(VEL < MINVEL) = MINVEL(VEL < MINVEL);
VEL(POS > MAXLIM) = (-1) .* VEL(POS > MAXLIM); % 反方向飞
POS(POS > MAXLIM) = MAXLIM(POS > MAXLIM);
VEL(POS < MINLIM) = (-1) .* VEL(POS < MINLIM); % 反方向飞
POS(POS < MINLIM) = MINLIM(POS < MINLIM);
end
function dom_vector = checkDomination(fitness)
Np = size(fitness, 1);
dom_vector = zeros(Np, 1);
all_perm = nchoosek(1 : Np, 2); % 1到Np中两两组合
all_perm = [all_perm; [all_perm(:, 2), all_perm(:, 1)]];
d = dominates(fitness(all_perm(:, 1), :), fitness(all_perm(:, 2), :));
dominated_particles = unique(all_perm(d == 1, 2));
dom_vector(dominated_particles) = 1;
end
function d = dominates(x, y)
% 判断支配
d = all(x <= y, 2) & any(x < y, 2);
end
function REP = updateGrid(REP, ngrid)
% 更新网格
ndim = size(REP.pos_fit, 2);
REP.hypercube_limits = zeros(ngrid + 1, ndim);
for dim = 1 : ndim
REP.hypercube_limits(:, dim) = linspace(min(REP.pos_fit(:, dim)), max(REP.pos_fit(:, dim)), ngrid + 1)';
end
% 计算优秀个体在网格中的分布个数
npar = size(REP.pos_fit, 1);
REP.grid_idx = zeros(npar, 1);
REP.grid_subidx = zeros(npar, ndim);
for n = 1 : 1 : npar
idnames = [];
for d = 1 : 1 : ndim
REP.grid_subidx(n, d) = find(REP.pos_fit(n, d) <= REP.hypercube_limits(:, d)', 1, 'first') - 1;
if(REP.grid_subidx(n,d) == 0)
REP.grid_subidx(n,d) = 1;
end
idnames = [idnames ',' num2str(REP.grid_subidx(n, d))];
end
REP.grid_idx(n) = eval(['sub2ind(ngrid.*ones(1,ndim)' idnames ');']);
end
REP.quality = zeros(ngrid, 2);
ids = unique(REP.grid_idx);
for i = 1 : length(ids)
REP.quality(i, 1) = ids(i); % 网格中个体个数越多质量越差
REP.quality(i, 2) = 10 / sum(REP.grid_idx == ids(i)); % 计算网格质量
end
end
function selected = selectLeader(REP)
% 选择网格中的一个个体做为最优个体
prob = cumsum(REP.quality(:,2));
sel_hyp = REP.quality(find(rand(1, 1) * max(prob) <= prob, 1, 'first'), 1);
idx = 1 : 1 : length(REP.grid_idx);
selected = idx(REP.grid_idx == sel_hyp);
selected = selected(randi(length(selected)));
end
function REP = deleteFromRepository(REP, n_extra, ngrid)
crowding = zeros(size(REP.pos, 1), 1);
for m = 1 : 1 : size(REP.pos_fit, 2)
[m_fit, idx] = sort(REP.pos_fit(:, m), 'ascend'); % 升序排列
m_up = [m_fit(2 : end); Inf];
m_down = [Inf; m_fit(1 : end - 1)];
distance = (m_up - m_down) ./ (max(m_fit) - min(m_fit));
[~, idx] = sort(idx, 'ascend');
crowding = crowding + distance(idx); % 计算距离
end
crowding(isnan(crowding)) = Inf;
[~, del_idx] = sort(crowding, 'ascend');
del_idx = del_idx(1 : n_extra);
REP.pos(del_idx, :) = [];
REP.pos_fit(del_idx, :) = [];
REP = updateGrid(REP, ngrid);
end
function POS = mutation(POS, gen, maxgen, Np, var_max, var_min, nVar, u_mut)
% 变异，避免陷入局部最优
% 个体分为3份
fract = Np / 3 - floor(Np / 3);
if fract < 0.5
sub_sizes =[ceil(Np / 3), round(Np / 3), round(Np / 3)];
else
sub_sizes =[round(Np / 3), round(Np / 3), floor(Np / 3)];
end
cum_sizes = cumsum(sub_sizes);
% 第二份变异
nmut = round(u_mut * sub_sizes(2));
if(nmut > 0)
idx = cum_sizes(1) + randperm(sub_sizes(2), nmut);
POS(idx, :) = repmat((var_max - var_min)', nmut, 1) .* rand(nmut, nVar) + repmat(var_min', nmut, 1);
end
% 第三份变异
per_mut = (1 - gen / maxgen) ^ (5 * nVar); % 变异率随代数变小
nmut = round(per_mut * sub_sizes(3));
if nmut > 0
idx = cum_sizes(2) + randperm(sub_sizes(3), nmut);
POS(idx, :) = repmat((var_max - var_min)', nmut, 1) .* rand(nmut, nVar) + repmat(var_min', nmut, 1);
end
end