Matlab遗传算法与非线性规划问题的结合
目录
准备
在编写自己的代码之前应该了解的是算法的基本结构,具体可参考我之前的文章大纲或者这里简单函数优化
文章虽然简单,但是也能了解基本该设置的参数与具体的算法结构。
遗传算法+非线性规划问题
在这篇文章以及之后的绝大多数文章里,遗传算法与各类算法的结合是重点,遗传算法自身的特点不再过多描述。遗传算法为其他算法进行服务,这也是在解决复杂问题时遗传算法的能力不足所导致的结果。在非线性规划中往往因为计算复杂,需要大量的时间去迭代计算,这个时间有时是一个工程无法承受的,因此采用遗传算法的提供较优解,框定解的起始点,减少非线性规划的计算时间,提高整体算法的计算速度。
更改了简单函数优化的目标,设定了此次的优化目标:
F
=
sin
(
10
π
⋅
x
y
)
+
sin
(
y
2
)
+
sin
(
5
z
)
F=\sin \left( 10\pi \cdot xy \right) +\sin \left( y^2 \right) +\sin \left( 5z \right)
F=sin(10π⋅xy)+sin(y2)+sin(5z)
在
x
,
y
,
z
∈
[
0
,
3
]
x,y,z\in \left[ 0,3 \right]
x,y,z∈[0,3]范围内求函数的最小值。
不难观察的是,函数值取得最大值3,
x
,
y
,
z
x,y,z
x,y,z需满足下列条件
{
10
π
⋅
x
y
=
k
1
π
2
⇒
x
y
=
k
1
20
y
2
=
k
2
π
2
⇒
y
=
k
2
π
2
5
z
=
k
3
π
2
⇒
z
=
k
3
π
10
x
,
y
,
z
∈
[
0
,
3
]
k
1
,
k
2
,
k
3
为
任
意
整
数
\left\{ \begin{array}{l} 10\pi \cdot xy=\dfrac{k_1\pi}{2}\Rightarrow xy=\dfrac{k_1}{20}\\ \\ y^2=\dfrac{k_2\pi}{2}\Rightarrow y=\sqrt{\dfrac{k_2\pi}{2}}\\ \\ 5z=\dfrac{k_3\pi}{2}\Rightarrow z=\dfrac{k_3\pi}{10}\\ \\ x,y,z∈\left[ 0,3 \right]\\ \end{array} \right. \,\,\,\,\,\,\,\,k_1,k_2,k_3为任意整数
⎩⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎧10π⋅xy=2k1π⇒xy=20k1y2=2k2π⇒y=2k2π5z=2k3π⇒z=10k3πx,y,z∈[0,3]k1,k2,k3为任意整数
m a t l a b matlab matlab自带的非线性规划函数 f m i n c o n fmincon fmincon
两者结合之前,
m
a
t
l
a
b
matlab
matlab拥有自带的非线性规划函数fmincon()
B
e
s
t
_
x
=
f
m
i
n
c
o
n
(
f
u
n
,
x
0
,
A
,
b
,
A
e
q
,
b
e
q
,
l
b
,
u
b
,
n
o
n
l
c
o
n
,
o
p
t
i
o
n
s
)
Best \_ x = fmincon(fun,x0,A,b,Aeq,beq,lb,ub,nonlcon,options)
Best_x=fmincon(fun,x0,A,b,Aeq,beq,lb,ub,nonlcon,options)
在这里"借鉴"
m
a
t
l
a
b
matlab
matlab中的解释
min
x
f
(
x
)
满足约束
{
c
(
x
)
≤
0
c
e
q
(
x
)
=
0
A
⋅
x
≤
b
A
e
q
⋅
x
=
b
e
q
l
b
≤
x
≤
u
b
\underset{x}{\min}\ f\left( x \right) \text{满足约束}\left\{ \begin{array}{l} c\left( x \right) \le 0\\ ceq\left( x \right) =0\\ A\cdot x\le b\\ Aeq\cdot x=beq\\ lb\le x\le ub\\ \end{array} \right.
xmin f(x)满足约束⎩⎪⎪⎪⎪⎨⎪⎪⎪⎪⎧c(x)≤0ceq(x)=0A⋅x≤bAeq⋅x=beqlb≤x≤ub
f
u
n
fun
fun为目标函数,求最小值,
x
0
x0
x0为起始点(也是遗传算法的目标)。
其中
c
(
x
)
,
c
e
q
(
x
)
c(x),ceq(x)
c(x),ceq(x)可以为非线性函数,例如简单的二次函数:
c
(
x
)
=
x
2
−
2
x
+
1
c(x)=x^2-2x+1
c(x)=x2−2x+1,两者均返回向量。
n
o
n
l
c
o
n
nonlcon
nonlcon参数同时包囊两者,常采用函数句柄输入
f
m
i
n
c
o
n
fmincon
fmincon。
其中
A
,
A
e
q
A,Aeq
A,Aeq为线性约束的自变量系数矩阵。
b
,
b
e
q
b,beq
b,beq同样为矩阵。
u
b
,
l
b
ub,lb
ub,lb为变量上下界。
值得注意的是,在约束条件里的符号是不能更改的。如果你的约束式符号与上述相反则需要不等式两边同时乘以-1
来改变符号。
由于本文中关于交叉与变异的操作方法与适应度无关,仅在选择操作与适应度有关,因此选择操作决定了遗传算法整体的选取值趋向(最大值还是最小值),且需要与非线性规划的趋向进行统一。
非线性寻优
- 采用 f m i n c o n fmincon fmincon进行非线性寻优
function ret = nonlinear(chrom,sizepop)
for i=1:sizepop
x=fmincon(@Fit_fun_nonliner,chrom(i,:)',[],[],[],[],[0 0 0 0 0],[5 5 5 5 5]);
ret(i,:)=x';
end
值得再强调的是,这里的非线性规划求得的是目标函数最小值。
实际上这里除了变量范围的限定并为对其他进行限制,具体的使用这里并非重点,有时间会特别出一期进行记录。
遗传算法
选择了谢菲尔德大学的Matlab遗传算法工具箱实现算法。
工具箱的导入极为简单,可以自行在网上下载导入。
由于与非线性规划结合,遗传算法在基本的种族的初始化,交叉,变异操作之后要对每个个体进行可行性检验,一般采用
w
h
i
l
e
while
while循环,直到操作产生的新个体符合要求。
可行性检验
- 对每一个个体进行范围检验
function flag=test(bound,code)
flag=1;
[n,m]=size(code);
for i=1:n
if code(i)<bound(i,1) || code(i)>bound(i,2) %范围检验
flag=0;
end
end
实际上,这里的约束不能过多,过多可能会导致程序困死于外层
w
h
i
l
e
while
while循环,这也是结合算法的缺点,如果约束过于复杂,遗传算法的交叉与变异可能无法产生合适的个体。
如果个体不符合范围要求则返回0,符合则返回1
种群初始化和计算初始种群适应度
- 初始参数设置
%% 清屏
clc
clear
close all
%% 定义遗传算法参数
MAXGEN=100; %进化代数
NIND=30; %种群规模
px=0.6; %交叉概率
pm=0.01; %变异概率
lenchrom=[1 1 1]; %变量字串长度
bound=[0 3;0 3;0 3]; %变量范围
- 设置种群结构体
- 产生满足一定条件的个体的种群
- 记录初始总体的最佳染色体
Chrom=struct('fitness',zeros(1,NIND), 'chrom',[]);
%随机初始化种群
for i=1:sizepop
Chrom.chrom(i,:)=Code(lenchrom,bound); %随机产生个体
x=Chrom.chrom(i,:);
Chrom.fitness(i)=Fit_fun_nonliner(x); %个体适应度
end
%找最好的染色体
[bestfitness,bestindex]=min(Chrom.fitness);
bestchrom=Chrom.chrom(bestindex,:); %最好的染色体
avgfitness=sum(Chrom.fitness)/NIND; %染色体的平均适应度
% 记录每一代进化中最好的适应度和平均适应度
trace=zeros(MAXGEN,2);
C o d e Code Code采用随机线性插值限制染色体在变量范围之内。
- C o d e Code Code函数展示如何添加染色体检验
function ret=Code(lenchrom,bound)
flag=0;
%检验不通过则重新进行插值
while flag==0
pick=rand(1,length(lenchrom)); %生成随机序列
ret=bound(:,1)'+(bound(:,2)-bound(:,1))'.*pick; %线性插值
flag=test(lenchrom,bound,ret); %检验染色体的可行性
end
实际上之后的交叉和变异操作与简单函数优化的区别并不大,仅仅添加了检验,不多做解释。
计算适应度
- 创建独立的.m文件计算适应度
function y = Fit_fun_nonliner(x)
y=-(sin(10*pi*x(1)*x(2))+sin(x(2)^2)+sin(5*x(3)));
end
这里稍微修改了简单函数优化的目标函数,添加负号求最大值。
个体选择
- 采用轮盘方法从种群中选择个体
Chrom=Select(Chrom,NIND); %选择
采用轮盘赌局的方法对种群进行选择,原轮盘选择法的选取趋向为最大值,因此在 S e l e c t Select Select函数中先对适应度取了倒数,适应度小的个体更容易被选取,与非线性规划保持统一。
交叉
- 采用单点交叉算子对个体的染色体进行重组
Chrom.chrom=Cross(px,lenchrom,Chrom.chrom,NIND,bound); %单点交叉
变异
- 对个体进行变异
Chrom.chrom=Mutation(pm,lenchrom,Chrom.chrom,NIND,[i MAXGEN],bound); %变异
非线性寻优
- 将当前进化的结果作为非线性规划的起始点进行优化
if mod(i,10)==0
Chrom.chrom=nonlinear(Chrom.chrom,NIND); %非线性寻优
end
每十代进行一次非线性寻优,求的在此起始点下的最优解
计算当前种群适应度记录最优解
- 计算当前所有个体的适应度。
- 找到并记录当代最优个体,更新最优解与平均解。
% 计算适应度
for j=1:NIND
x=Chrom.chrom(j,:);
Chrom.fitness(j)=Fit_fun_nonliner(x);
end
%找到最小和最大适应度的染色体及它们在种群中的位置
[newbestfitness,newbestindex]=min(Chrom.fitness);
[worestfitness,worestindex]=max(Chrom.fitness);
% 代替上一次进化中最好的染色体
if bestfitness>newbestfitness
bestfitness=newbestfitness;
bestchrom=Chrom.chrom(newbestindex,:);
end
Chrom.chrom(worestindex,:)=bestchrom;
Chrom.fitness(worestindex)=bestfitness;
avgfitness=sum(Chrom.fitness)/NIND;
trace(i,:)=[avgfitness,bestfitness]; %记录每一代进化中最好的适应度和平均适应度
当平均解与最优解重合时说明当前种群中大多均为最优解,在一定程度代表当前种群已经较为稳定
算法计算结果
虽然在40代左右种群有不稳定的现象,但随着迭代次数再次增加,最终种群中绝大多数个体均接近最优个体。
说明现在的种群已经足够适应当前的"环境",此时已经产生了较好解,最优值与最优解如下:
最优值 | x x x | y y y | z z z |
---|---|---|---|
3.0000 3.0000 3.0000 | 0.0399 0.0399 0.0399 | 1.2533 1.2533 1.2533 | 0.3142 0.3142 0.3142 |
对结果进行验证与预测结果一致。
x y xy xy | y 2 y^2 y2 | 5 z 5z 5z |
---|---|---|
0.050 0.050 0.050 | 1.571 1.571 1.571 | 1.571 1.571 1.571 |
1 / 2 1/2 1/2 | π / 2 \pi/2 π/2 | π / 2 \pi/2 π/2 |
在这里只是举了一个简单的例子,哪怕不使用非线性规划,单采用遗传算法也能有较好的结果。重点在于两者的结合,具体问题要具体分析。
代码地址
完整代码在 GitHub网站上,如果有条件可以支持一下。现在关注我你就是老粉了
参考文章
本文内容参考下列文章,本篇使用的函数与参考文章基本一致,只做了小幅度更改。如有任何问题还请联系该账号
史峰. MATLAB智能算法30个案例分析[M]. 北京航空航天大学出版社, 2011.