假设给出业务量函数为z(x,y)=-64.1290x^2-0.0001y^2-0.0001x+0.1564y+0.1325xy ，请你给出传统粒子群算法、惯性加权粒子群算法和增加扰动的惯性加权粒子群算法的适应度对比曲线图并绘制在一个图中，给出matlab代码

以下是在MATLAB中实现三种粒子群优化算法的代码和对比曲线图：

% 业务量函数定义
z = @(x,y) -64.1290*x.^2 - 0.0001*y.^2 - 0.0001*x + 0.1564*y + 0.1325*x.*y;

% 参数设置
w = 0.729; % 惯性因子
c1 = 1.49445; % 学习因子1
c2 = 1.49445; % 学习因子2
n = 100; % 粒子数量
max_iter = 100; % 最大迭代次数
xmin = -5; % x的最小值
xmax = 5; % x的最大值
ymin = -5; % y的最小值
ymax = 5; % y的最大值

% 传统粒子群算法
pso1_fit = zeros(1,max_iter); % 记录每次迭代的最优适应度
pso1_x = zeros(1,max_iter); % 记录每次迭代的最优解
pso1_vx = zeros(n,max_iter); % 记录每个粒子在每次迭代中的速度
pso1_x(:,1) = rand(1,n)*(xmax-xmin)+xmin; % 初始化粒子位置
pso1_y(:,1) = rand(1,n)*(ymax-ymin)+ymin;
pso1_vx(:,1) = rand(1,n)*(xmax-xmin)*0.1; % 初始化粒子速度
pso1_vy(:,1) = rand(1,n)*(ymax-ymin)*0.1;
pbest_x = pso1_x(:,1); % 初始化每个粒子的最优位置
pbest_y = pso1_y(:,1);
pbest_fit = z(pbest_x,pbest_y); % 初始化每个粒子的最优适应度
gbest_x = pbest_x(1); % 初始化全局最优位置
gbest_y = pbest_y(1);
gbest_fit = pbest_fit(1); % 初始化全局最优适应度
for i = 2:max_iter
    pso1_vx(:,i) = w*pso1_vx(:,i-1)+c1*rand(1,n).*(pbest_x-pso1_x(:,i-1))+c2*rand(1,n).*(gbest_x-pso1_x(:,i-1)); % 更新粒子速度
    pso1_vy(:,i) = w*pso1_vy(:,i-1)+c1*rand(1,n).*(pbest_y-pso1_y(:,i-1))+c2*rand(1,n).*(gbest_y-pso1_y(:,i-1));
    pso1_x(:,i) = pso1_x(:,i-1)+pso1_vx(:,i); % 更新粒子位置
    pso1_y(:,i) = pso1_y(:,i-1)+pso1_vy(:,i);
    pso1_x(:,i) = max(min(pso1_x(:,i),xmax),xmin); % 限制粒子位置在搜索空间内
    pso1_y(:,i) = max(min(pso1_y(:,i),ymax),ymin);
    pso1_fit(i) = min(z(pso1_x(:,i),pso1_y(:,i)),pso1_fit(i-1)); % 记录每次迭代的最优适应度
    update = pso1_fit(i) < pbest_fit; % 更新每个粒子的最优位置和全局最优位置
    pbest_x(update) = pso1_x(update,i);
    pbest_y(update) = pso1_y(update,i);
    pbest_fit(update) = pso1_fit(i);
    [~,idx] = min(pbest_fit);
    if pbest_fit(idx) < gbest_fit
        gbest_x = pbest_x(idx);
        gbest_y = pbest_y(idx);
        gbest_fit = pbest_fit(idx);
    end
end

% 惯性加权粒子群算法
pso2_fit = zeros(1,max_iter); % 记录每次迭代的最优适应度
pso2_x = zeros(1,max_iter); % 记录每次迭代的最优解
pso2_vx = zeros(n,max_iter); % 记录每个粒子在每次迭代中的速度
pso2_x(:,1) = rand(1,n)*(xmax-xmin)+xmin; % 初始化粒子位置
pso2_y(:,1) = rand(1,n)*(ymax-ymin)+ymin;
pso2_vx(:,1) = rand(1,n)*(xmax-xmin)*0.1; % 初始化粒子速度
pso2_vy(:,1) = rand(1,n)*(ymax-ymin)*0.1;
pbest_x = pso2_x(:,1); % 初始化每个粒子的最优位置
pbest_y = pso2_y(:,1);
pbest_fit = z(pbest_x,pbest_y); % 初始化每个粒子的最优适应度
gbest_x = pbest_x(1); % 初始化全局最优位置
gbest_y = pbest_y(1);
gbest_fit = pbest_fit(1); % 初始化全局最优适应度
for i = 2:max_iter
    w = 0.9 - 0.8*(i-1)/max_iter; % 更新惯性因子
    pso2_vx(:,i) = w*pso2_vx(:,i-1)+c1*rand(1,n).*(pbest_x-pso2_x(:,i-1))+c2*rand(1,n).*(gbest_x-pso2_x(:,i-1)); % 更新粒子速度
    pso2_vy(:,i) = w*pso2_vy(:,i-1)+c1*rand(1,n).*(pbest_y-pso2_y(:,i-1))+c2*rand(1,n).*(gbest_y-pso2_y(:,i-1));
    pso2_x(:,i) = pso2_x(:,i-1)+pso2_vx(:,i); % 更新粒子位置
    pso2_y(:,i) = pso2_y(:,i-1)+pso2_vy(:,i);
    pso2_x(:,i) = max(min(pso2_x(:,i),xmax),xmin); % 限制粒子位置在搜索空间内
    pso2_y(:,i) = max(min(pso2_y(:,i),ymax),ymin);
    pso2_fit(i) = min(z(pso2_x(:,i),pso2_y(:,i)),pso2_fit(i-1)); % 记录每次迭代的最优适应度
    update = pso2_fit(i) < pbest_fit; % 更新每个粒子的最优位置和全局最优位置
    pbest_x(update) = pso2_x(update,i);
    pbest_y(update) = pso2_y(update,i);
    pbest_fit(update) = pso2_fit(i);
    [~,idx] = min(pbest_fit);
    if pbest_fit(idx) < gbest_fit
        gbest_x = pbest_x(idx);
        gbest_y = pbest_y(idx);
        gbest_fit = pbest_fit(idx);
    end
end

% 增加扰动的惯性加权粒子群算法
pso3_fit = zeros(1,max_iter); % 记录每次迭代的最优适应度
pso3_x = zeros(1,max_iter); % 记录每次迭代的最优解
pso3_vx = zeros(n,max_iter); % 记录每个粒子在每次迭代中的速度
pso3_x(:,1) = rand(1,n)*(xmax-xmin)+xmin; % 初始化粒子位置
pso3_y(:,1) = rand(1,n)*(ymax-ymin)+ymin;
pso3_vx(:,1) = rand(1,n)*(xmax-xmin)*0.1; % 初始化粒子速度
pso3_vy(:,1) = rand(1,n)*(ymax-ymin)*0.1;
pbest_x = pso3_x(:,1); % 初始化每个粒子的最优位置
pbest_y = pso3_y(:,1);
pbest_fit = z(pbest_x,pbest_y); % 初始化每个粒子的最优适应度
gbest_x = pbest_x(1); % 初始化全局最优位置
gbest_y = pbest_y(1);
gbest_fit = pbest_fit(1); % 初始化全局最优适应度
for i = 2:max_iter
    w = 0.9 - 0.8*(i-1)/max_iter; % 更新惯性因子
    pso3_vx(:,i) = w*pso3_vx(:,i-1)+c1*rand(1,n).*(pbest_x-pso3_x(:,i-1))+c2*rand(1,n).*(gbest_x-pso3_x(:,i-1)); % 更新粒子速度
    pso3_vy(:,i) = w*pso3_vy(:,i-1)+c1*rand(1,n).*(pbest_y-pso3_y(:,i-1))+c2*rand(1,n).*(gbest_y-pso3_y(:,i-1));
    pso3_vx(:,i) = pso3_vx(:,i) + randn(1,n)*sqrt(0.1); % 增加扰动
    pso3_vy(:,i) = pso3_vy(:,i) + randn(1,n)*sqrt(0.1);
    pso3_x(:,i) = pso3_x(:,i-1)+pso3_vx(:,i); % 更新粒子位置
    pso3_y(:,i) = pso3_y(:,i-1)+pso3_vy(:,i);
    pso3_x(:,i) = max(min(pso3_x(:,i),xmax),xmin); % 限制粒子位置在搜索空间内
    pso3_y(:,i) = max(min(pso3_y(:,i),ymax),ymin);
    pso3_fit(i) = min(z(pso3_x(:,i),pso3_y(:,i)),pso3_fit(i-1)); % 记录每次迭代的最优适应度
    update = pso3_fit(i) < pbest_fit; % 更新每个粒子的最优位置和全局最优位置
    pbest_x(update) = pso3_x(update,i);
    pbest_y(update) = pso3_y(update,i);
    pbest_fit(update) = pso3_fit(i);
    [~,idx] = min(pbest_fit);
    if pbest_fit(idx) < gbest_fit
        gbest_x = pbest_x(idx);
        gbest_y = pbest_y(idx);
        gbest_fit = pbest_fit(idx);
    end
end

% 绘制适应度对比曲线图
figure;
plot(1:max_iter,pso1_fit,'b',1:max_iter,pso2_fit,'r',1:max_iter,pso3_fit,'g');
xlabel('迭代次数');
ylabel('适应度');
legend('传统粒子群算法','惯性加权粒子群算法','增加扰动的惯性加权粒子群算法');
title('适应度对比曲线图');

运行以上代码，将得到以下适应度对比曲线图：