gafit2.m

malab语言实现bspline曲线拟合

clear all;
close all;
clc;

NIND = 40;
MAXGEN = 200;
m = 34;
p = 3;
W = .9;
C1 = 2;
C2 = 2;
Max_velocity = 0.02;
Max_Nv = 13;

Q = [  
  0.000000  0.000000
  0.012500  0.019300
  0.025000  0.031700
  0.050000  0.051300
  0.075000  0.066400
  0.100000  0.078000
  0.150000  0.093400
  0.200000  0.101300
  0.250000  0.104400
  0.300000  0.104800
  0.400000  0.100200
  0.500000  0.090500
  0.600000  0.077100
  0.700000  0.061000
  0.800000  0.042800
  0.900000  0.022900
  0.950000  0.012400
  1.000000  0.001600
  1.000000 -0.001600
  0.950000  0.004400
  0.900000  0.010000
  0.800000  0.020200
  0.700000  0.028600
  0.600000  0.035000
  0.500000  0.038600
  0.400000  0.038500
  0.300000  0.033300
  0.250000  0.028200
  0.200000  0.021700
  0.150000  0.014500
  0.100000  0.006800
  0.075000  0.002800
  0.050000 -0.001000
  0.025000 -0.004200
  0.012500 -0.005000];
 
[chrom,v,pbest_chrom,pbest,gbest_chrom,gbest] = initialpop2(NIND,m,p,Q);    
[objv1,fitnv1] = objfun2(m,p,chrom,Q);
gen = 1;

while gen < MAXGEN   

    chrom = update_position(chrom,v,NIND,m);

    % ================寻找pbest和gbest==========    
    [objv,fitnv] = objfun2(m,p,chrom,Q);
    for i = 1:NIND
        if fitnv(i) > pbest(i)
            pbest(i) = fitnv(i);
            pbest_chrom(i,:) = chrom(i,:);
        end
    end
    [fitnvmax,index_gbest] = max(fitnv);    

    if fitnvmax > gbest
        gbest = fitnvmax;
        gbest_chrom = chrom(index_gbest,:);
        draw(m,p,gbest_chrom,Q);
    end
    
    % ===========================================

    v = update_velocity(NIND,chrom,W,C1,C2,Max_velocity,Max_Nv,pbest_chrom,gbest_chrom,v,m);

    trace(1,gen) = fitnvmax;
    trace(2,gen) = sum(fitnv)/size(fitnv,1);
    trace(3,gen)=chrom(index_gbest,m+2);
    [ssemin,I]=min(objv);
    trace(4,gen)=chrom(I,m+2);
    gen = gen + 1

    if min(objv) < 0.0001
        gen
        gen = MAXGEN;
    end
    pause(0.01)   
end


% ============最后的输出==================
% disp(['RMSE0 = ',num2str(min(objv1))]);
% disp(['Fitnv0 = ',num2str(max(fitnv1))]);
% disp(['RMSE1 = ',num2str(min(objv))]);
% disp(['Fitnv1 = ',num2str(max(fitnv))]);

if((min(objv1)>min(objv))&&(max(fitnv1)<max(fitnv)))
    disp('Successful !');
    disp(['best N = ',num2str(gbest_chrom(m+2))]);
%   disp(['N = ',num2str(chrom(:,m+2)')]);
    [a,b] = size(find(chrom(:,m+2) == gbest_chrom(m+2)));
    disp(['number of best_N = ',num2str(a)]);
    
    figure; 
    plot(trace(1,:)');
    hold on;
    plot(trace(2,:)','-.');
    grid; 
    figure; 
    plot(trace(3,:)');
    hold on;
    plot(trace(4,:)','r-');
    xlabel('generation');
    ylabel('control point number');
    legend('适应度最大曲线的控制顶点数','误差平方和最小曲线的控制顶点数')
    grid; 
else
    disp('Failed !');
end