目录
前言
对于非线性回归拟合曲线，matlab中有函数fmincon，但是此函数适用于求解局部最优，难以找到全局最优，所以需要结合遗传算法或模拟退火等算法来寻求全局最优解。
1、非线性寻优
非线性寻优需要确定目标函数，约束条件，初值，上下限，利用fmincon函数直接求解。
0.0 原始函数
%% 含未知数函数
function y = fun33(x,A,B,C,D)
y = A*sin(x)+B*x.^2+C*x+D*log(x);
end1.1 目标函数
取19个点进行曲线拟合，目标函数为拟合曲线与散点的差值的绝对值。
%% 目标函数
function f=fun11(x)
A = x(1);
B = x(2);
C = x(3);
D = x(4);
x_0 = 1:0.5:10;
y_0 = [5.84147098480790 11.6193554190367 17.6818861490655 24.2636350716006 31.5355691627323 39.6602686462919 48.7883749491716 59.0387794694400 70.4788273750733 83.1134520433833 96.8876223787133 111.702328695694 127.440627194940 143.997612058944 161.307124413343 179.358751766609 198.201016794587 217.930016073964 238.666319261087];
f = 0;
for i = 1:length(x_0)
y = A*sin(x_0(i))+B*x_0(i).^2+C*x_0(i)+D*log(x_0(i));
f = f + abs(y-y_0(i));
end
end1.2 约束条件
有等式约束与不等式约束，如果是线性约束直接写到fmincon函数中有A、Aeq、B、Beq输入参数，但是有很多情况不是线性约束而是非线性约束，如谐振时LC=1/(w)^2，所以要建立约束条件函数。
%% 约束条件
function [g,h]=fun22(x)
g=[];%不等式约束<=0
h=[x(1)*x(2)-2];%等式约束
end1.3 初值、上下限设置
%% 计算
%初始化
options = optimoptions('fmincon','MaxFunEvals',1000000,'TolFun',10^(-100000),'Algorithm','sqp');
%设定初值（非常关键，决定了结果的准确性与计算的速度）
x0 = [0.5 0.5 0.5 0.5];1.4 调用fimincon函数并计算
%得到x y
[x,y] = fmincon(@fun11,x0,[],[],[],[],[0 0 0 0],Inf,@fun22);1.5 结果显示
A = x(1);
B = x(2);
C = x(3);
D = x(4);
%求解得到的函数
x_0 = 1:0.5:10;
yy = fun33(x_0,A,B,C,D);
%原始数据
y_0 = [5.84147098480790 11.6193554190367 17.6818861490655 24.2636350716006 31.5355691627323 39.6602686462919 48.7883749491716 59.0387794694400 70.4788273750733 83.1134520433833 96.8876223787133 111.702328695694 127.440627194940 143.997612058944 161.307124413343 179.358751766609 198.201016794587 217.930016073964 238.666319261087];
%画图
figure;
plot(x_0,y_0,'--','LineWidth',2);
hold on;
plot(x_0,yy,'r-');
legend('RAW_{data}','FIT_{data}','FontSize',25);%latax?对结果进行显示，最终效果如下图，误差为1e-4，初步设置ABCD为1，2，3，4，最终拟合结果ABCD为1，2，3，3.999899107223816。
2、遗传算法+非线性寻优
既然fimincon函数这么准了，那为什么还有使用遗传算法进行寻优呢？因为fimincon函数有局限性，此函数在解决复杂非线性问题时容易陷入局部最优解，所以需要使用遗传算法、模拟退火算法等来寻找全局最优解。
2.1 建立目标函数、约束条件、原始函数
%% 目标函数
function f=fun11(x)
A = x(1);
B = x(2);
C = x(3);
D = x(4);
x_0 = 1:0.5:10;
y_0 = [5.84147098480790 11.6193554190367 17.6818861490655 24.2636350716006 31.5355691627323 39.6602686462919 48.7883749491716 59.0387794694400 70.4788273750733 83.1134520433833 96.8876223787133 111.702328695694 127.440627194940 143.997612058944 161.307124413343 179.358751766609 198.201016794587 217.930016073964 238.666319261087];
f = 0;
for i = 1:length(x_0)
y = A*sin(x_0(i))+B*x_0(i).^2+C*x_0(i)+D*log(x_0(i));
f = f + abs(y-y_0(i));
end
end
%% 约束条件
function [g,h]=fun22(x)
g=[];%不等式约束<=0
h=[x(1)*x(2)-2];%等式约束
end
%% 含未知数函数
function y = fun33(x,A,B,C,D)
y = A*sin(x)+B*x.^2+C*x+D*log(x);
end2.2 遗传算法初始化
对遗传算法的一些参数进行初始化。
%% 遗传算法参数
maxgen=30; %进化代数
sizepop=100; %种群规模
pcross=[0.6]; %交叉概率
pmutation=[0.01]; %变异概率
lenchrom=[1 1 1 1]; %变量字串长度
bound=[0 15;0 15;0 15;0 15]; %变量范围2.3 个体初始化
初始化个体计算群体
%% 个体初始化
individuals=struct('fitness',zeros(1,sizepop), 'chrom',[]); %种群结构体
avgfitness=[]; %种群平均适应度
bestfitness=[]; %种群最佳适应度
bestchrom=[]; %适应度最好染色体
% 初始化种群
for i=1:sizepop
individuals.chrom(i,:)=Code(lenchrom,bound); %随机产生个体
x=individuals.chrom(i,:);
individuals.fitness(i)=fun11(x); %个体适应度
end
%找最好的染色体
[bestfitness, bestindex]=min(individuals.fitness);
bestchrom=individuals.chrom(bestindex,:); %最好的染色体
avgfitness=sum(individuals.fitness)/sizepop; %染色体的平均适应度
% 记录每一代进化中最好的适应度和平均适应度
trace=[];2.4 进行进化
%% 进化开始
for i=1:maxgen
% 选择操作
individuals=Select(individuals,sizepop);
avgfitness=sum(individuals.fitness)/sizepop;
% 交叉操作
individuals.chrom=Cross(pcross,lenchrom,individuals.chrom,sizepop,bound);
% 变异操作
individuals.chrom=Mutation(pmutation,lenchrom,individuals.chrom,sizepop,[i maxgen],bound);
%每进化10代，以所得值为初始值进行非线性寻优
if mod(i,10)==0
individuals.chrom=nonlinear(individuals.chrom,sizepop);
end
% 计算适应度
for j=1:sizepop
x=individuals.chrom(j,:);
individuals.fitness(j)=fun11(x);
end
%找到最小和最大适应度的染色体及它们在种群中的位置
[newbestfitness,newbestindex]=min(individuals.fitness);
[worestfitness,worestindex]=max(individuals.fitness);
% 代替上一次进化中最好的染色体
if bestfitness>newbestfitness
bestfitness=newbestfitness;
bestchrom=individuals.chrom(newbestindex,:);
end
individuals.chrom(worestindex,:)=bestchrom;
individuals.fitness(worestindex)=bestfitness;
avgfitness=sum(individuals.fitness)/sizepop;
trace=[trace;avgfitness bestfitness]; %记录每一代进化中最好的适应度和平均适应度
end
%进化结束2.5 结果显示
显示出最佳函数值与x值，并且迭代过程可视化，最终拟合结果可视化。
%% 结果显示
figure
[r, c]=size(trace);
plot([1:r]',trace(:,1),'r-',[1:r]',trace(:,2),'b--');
title(['函数值曲线 ' '终止代数＝' num2str(maxgen)],'fontsize',12);
xlabel('进化代数','fontsize',12);ylabel('函数值','fontsize',12);
legend('各代平均值','各代最佳值','fontsize',12);
disp('函数值 变量');
% 窗口显示
disp([bestfitness x]);
grid on
% figure(2)
A = x(1);
B = x(2);
C = x(3);
D = x(4);
%求解得到的函数
x_0 = 1:0.5:10;
yy = fun33(x_0,A,B,C,D);
%原始数据
y_0 = [5.84147098480790 11.6193554190367 17.6818861490655 24.2636350716006 31.5355691627323 39.6602686462919 48.7883749491716 59.0387794694400 70.4788273750733 83.1134520433833 96.8876223787133 111.702328695694 127.440627194940 143.997612058944 161.307124413343 179.358751766609 198.201016794587 217.930016073964 238.666319261087];
%画图
figure;
plot(x_0,y_0,'--','LineWidth',2);
hold on;
plot(x_0,yy,'r-');
legend('RAW_{data}','FIT_{data}','FontSize',25);%latax?2.6 编码
function ret=Code(lenchrom,bound)
%本函数将变量编码成染色体，用于随机初始化一个种群
% lenchrom input : 染色体长度
% bound input : 变量的取值范围
% ret output: 染色体的编码值
flag=0;
while flag==0
pick=rand(1,length(lenchrom));
ret=bound(:,1)'+(bound(:,2)-bound(:,1))'.*pick; %线性插值
flag=test(lenchrom,bound,ret); %检验染色体的可行性
end
end2.7 编码检验约束判断
%% 判断约束
function flag=test(lenchrom,bound,code)
% lenchrom input : 染色体长度
% bound input : 变量的取值范围
% code output: 染色体的编码值
flag=1;
[n,m]=size(code);
for i=1:n
% if code(i)<bound(i,1) || code(i)>bound(i,2)||code(3)*code(4)<1/(2*330*pi)^2||code(3)*code(4)>1/(2*300*pi)^2
if code(i)<bound(i,1) || code(i)>bound(i,2)
flag=0;
end
end
end2.8 非线性寻优
%% 非线性寻优
function ret = nonlinear(chrom,sizepop)
for i=1:sizepop
x=fmincon(@fun11,chrom(i,:)',[],[],[],[],[0 0 0 0],Inf,@fun22);
ret(i,:)=x';
end
end2.9 变异
function ret=Mutation(pmutation,lenchrom,chrom,sizepop,pop,bound)
% 本函数完成变异操作
% pcorss input : 变异概率
% lenchrom input : 染色体长度
% chrom input : 染色体群
% sizepop input : 种群规模
% pop input : 当前种群的进化代数和最大的进化代数信息
% ret output : 变异后的染色体
for i=1:sizepop
% 随机选择一个染色体进行变异
pick=rand;
while pick==0
pick=rand;
end
index=ceil(pick*sizepop);
% 变异概率决定该轮循环是否进行变异
pick=rand;
if pick>pmutation
continue;
end
flag=0;
while flag==0
% 变异位置
pick=rand;
while pick==0
pick=rand;
end
pos=ceil(pick*sum(lenchrom)); %随机选择了染色体变异的位置，即选择了第pos个变量进行变异
v=chrom(i,pos);
v1=v-bound(pos,1);
v2=bound(pos,2)-v;
pick=rand; %变异开始
if pick>0.5
delta=v2*(1-pick^((1-pop(1)/pop(2))^2));
chrom(i,pos)=v+delta;
else
delta=v1*(1-pick^((1-pop(1)/pop(2))^2));
chrom(i,pos)=v-delta;
end %变异结束
flag=test(lenchrom,bound,chrom(i,:)); %检验染色体的可行性
end
end
ret=chrom;
end2.10 交叉
function ret=Cross(pcross,lenchrom,chrom,sizepop,bound)
%本函数完成交叉操作
% pcorss input : 交叉概率
% lenchrom input : 染色体的长度
% chrom input : 染色体群
% sizepop input : 种群规模
% ret output : 交叉后的染色体
for i=1:sizepop
% 随机选择两个染色体进行交叉
pick=rand(1,2);
while prod(pick)==0
pick=rand(1,2);
end
index=ceil(pick.*sizepop);
% 交叉概率决定是否进行交叉
pick=rand;
while pick==0
pick=rand;
end
if pick>pcross
continue;
end
flag=0;
while flag==0
% 随机选择交叉位置
pick=rand;
while pick==0
pick=rand;
end
pos=ceil(pick.*sum(lenchrom)); %随机选择进行交叉的位置，即选择第几个变量进行交叉，注意：两个染色体交叉的位置相同
pick=rand; %交叉开始
v1=chrom(index(1),pos);
v2=chrom(index(2),pos);
chrom(index(1),pos)=pick*v2+(1-pick)*v1;
chrom(index(2),pos)=pick*v1+(1-pick)*v2; %交叉结束
flag1=test(lenchrom,bound,chrom(index(1),:)); %检验染色体1的可行性
flag2=test(lenchrom,bound,chrom(index(2),:)); %检验染色体2的可行性
if flag1*flag2==0
flag=0;
else flag=1;
end %如果两个染色体不是都可行，则重新交叉
end
end
ret=chrom;
end2.11 选择
function ret=Select(individuals,sizepop)
% 本函数对每一代种群中的染色体进行选择，以进行后面的交叉和变异
% individuals input : 种群信息
% sizepop input : 种群规模
% opts input : 选择方法的选择
% ret output : 经过选择后的种群
individuals.fitness= 1./(individuals.fitness);
sumfitness=sum(individuals.fitness);
sumf=individuals.fitness./sumfitness;
index=[];
for i=1:sizepop %转sizepop次轮盘
pick=rand;
while pick==0
pick=rand;
end
for j=1:sizepop
pick=pick-sumf(j);
if pick<0
index=[index j];
break; %寻找落入的区间，此次转轮盘选中了染色体i，注意：在转sizepop次轮盘的过程中，有可能会重复选择某些染色体
end
end
end
individuals.chrom=individuals.chrom(index,:);
individuals.fitness=individuals.fitness(index);
ret=individuals;
end计算结果如下图：
看起来和非线性寻优没有差别？非也，得到的误差为1e-6，而得到的ABCD为1.000000000552358 1.999999997526807 3.000000158778851 3.999999450511331
看起来这四个数不太准，对吗？非也，因为原始数据是在ABCD为[1 2 3 4]时得到的，但是一定能保证得到的数值一定是准确的吗？不能，因为计算机有有精度，有精度就一定有误差，所以遗传算法得到的数值是基于“不准确”的最准确，如果最终得到ABCD为1234反而不准确。
3、参考资料
【MATLAB】求解约束条件下的目标函数最值（fmincon用法解析）_fmincon函数用法-CSDN博客
遗传算法与非线性规划的结合算法（附代码）_sumfitness-CSDN博客
2148
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
