MATLAB中 f1=@(x,y,z) x+y+z; if (x.*y.*z-(x+y+z).*(x.*y+y.*z+z.*x)./3+(2.*(x+y+z).^3)./27) >= 0 f21=@(x,y,z) ((((((x+y+z).^2))/3-(x.*y+y.*z+z.*x))).^(3/2)-... b1.*(x.*y.*z-(x+y+z).*(x.*y+y.*z+z.*x)./3+(2.*(x+y+z).^3)./27))-yield; else f22=@(x,y,z) ((((((x+y+z).^2))/3-(x.*y+y.*z+z.*x))).^(3/2)-... b2.*(x.*y.*z-(x+y+z).*(x.*y+y.*z+z.*x)./3+(2.*(x+y+z).^3)./27))-yield; end f2 = [f21 f22] 这句语法对吗

这个语法是不正确的，因为 f21 和 f22 是在 if-else 语句块中定义的局部函数句柄，不能在 if-else 语句块外部使用。如果需要在 if-else 语句块外部使用 f2，可以先在外部定义 f21 和 f22 函数句柄，然后在 if-else 语句块中分别给 f2 赋值，如下所示：

f21 = @(x,y,z) ((((((x+y+z).^2))/3-(x.*y+y.*z+z.*x))).^(3/2)-b1.*(x.*y.*z-(x+y+z).*(x.*y+y.*z+z.*x)./3+(2.*(x+y+z).^3)./27))-yield;
f22 = @(x,y,z) ((((((x+y+z).^2))/3-(x.*y+y.*z+z.*x))).^(3/2)-b2.*(x.*y.*z-(x+y+z).*(x.*y+y.*z+z.*x)./3+(2.*(x+y+z).^3)./27))-yield;
if (x.*y.*z-(x+y+z).*(x.*y+y.*z+z.*x)./3+(2.*(x+y+z).^3)./27) >= 0
    f2 = f21;
else
    f2 = f22;
end

这样就可以正确地定义 f2。

相关问题

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

传统粒子群算法：

% 传统PSO算法
clear;clc;
N = 50;  % 粒子群数目
maxgen = 100;   % 迭代次数
w = 0.8;    % 惯性权重
c1 = 2;     % 自我认知学习因子
c2 = 2;     % 社会认知学习因子
r1 = rand(N,1);
r2 = rand(N,1);
Vmax = 10;  % 最大速度
Vmin = -10; % 最小速度
% 初始化种群位置和速度
for i=1:N
    x(i) = rand(1)*100;
    y(i) = rand(1)*100;
    Vx(i) = (rand(1)-0.5)*2*10;
    Vy(i) = (rand(1)-0.5)*2*10;
end
% 计算适应度值
for i=1:N
    z(i) = -64.1290*x(i)^2-0.0001*y(i)^2-0.0001*x(i)+0.1564*y(i)+0.1325*x(i)*y(i);
end
% 初始化全局最优解
[bestval,idx] = min(z);
pbestx = x;
pbesty = y;
gbestx = x(idx);
gbesty = y(idx);
% 迭代寻优
for gen=1:maxgen
    % 更新速度和位置
    for i=1:N
        Vx(i) = w*Vx(i)+c1*r1(i)*(pbestx(i)-x(i))+c2*r2(i)*(gbestx-x(i));
        Vy(i) = w*Vy(i)+c1*r1(i)*(pbesty(i)-y(i))+c2*r2(i)*(gbesty-y(i));
        % 限制速度范围
        Vx(i) = min(max(Vx(i),Vmin),Vmax);
        Vy(i) = min(max(Vy(i),Vmin),Vmax);
        % 更新位置
        x(i) = x(i)+Vx(i);
        y(i) = y(i)+Vy(i);
        % 限制位置范围
        x(i) = min(max(x(i),0),100);
        y(i) = min(max(y(i),0),100);
    end
    % 计算适应度值
    for i=1:N
        z(i) = -64.1290*x(i)^2-0.0001*y(i)^2-0.0001*x(i)+0.1564*y(i)+0.1325*x(i)*y(i);
    end
    % 更新个体最优解和全局最优解
    for i=1:N
        if z(i) < -64.1290*pbestx(i)^2-0.0001*pbesty(i)^2-0.0001*pbestx(i)+0.1564*pbesty(i)+0.1325*pbestx(i)*pbesty(i)
            pbestx(i) = x(i);
            pbesty(i) = y(i);
        end
        if z(i) < bestval
            bestval = z(i);
            gbestx = x(i);
            gbesty = y(i);
        end
    end
    % 输出迭代结果
    fprintf('第%d代，最优解：%f，坐标（%f，%f）\n',gen,bestval,gbestx,gbesty);
end

惯性加权粒子群算法：

% 惯性加权PSO算法
clear;clc;
N = 50;  % 粒子群数目
maxgen = 100;   % 迭代次数
wmax = 0.9; % 最大惯性权重
wmin = 0.4; % 最小惯性权重
c1 = 2;     % 自我认知学习因子
c2 = 2;     % 社会认知学习因子
r1 = rand(N,1);
r2 = rand(N,1);
Vmax = 10;  % 最大速度
Vmin = -10; % 最小速度
% 初始化种群位置和速度
for i=1:N
    x(i) = rand(1)*100;
    y(i) = rand(1)*100;
    Vx(i) = (rand(1)-0.5)*2*10;
    Vy(i) = (rand(1)-0.5)*2*10;
end
% 计算适应度值
for i=1:N
    z(i) = -64.1290*x(i)^2-0.0001*y(i)^2-0.0001*x(i)+0.1564*y(i)+0.1325*x(i)*y(i);
end
% 初始化全局最优解
[bestval,idx] = min(z);
pbestx = x;
pbesty = y;
gbestx = x(idx);
gbesty = y(idx);
% 迭代寻优
for gen=1:maxgen
    % 计算当前惯性权重
    w = wmax - (wmax-wmin)*gen/maxgen;
    % 更新速度和位置
    for i=1:N
        Vx(i) = w*Vx(i)+c1*r1(i)*(pbestx(i)-x(i))+c2*r2(i)*(gbestx-x(i));
        Vy(i) = w*Vy(i)+c1*r1(i)*(pbesty(i)-y(i))+c2*r2(i)*(gbesty-y(i));
        % 限制速度范围
        Vx(i) = min(max(Vx(i),Vmin),Vmax);
        Vy(i) = min(max(Vy(i),Vmin),Vmax);
        % 更新位置
        x(i) = x(i)+Vx(i);
        y(i) = y(i)+Vy(i);
        % 限制位置范围
        x(i) = min(max(x(i),0),100);
        y(i) = min(max(y(i),0),100);
    end
    % 计算适应度值
    for i=1:N
        z(i) = -64.1290*x(i)^2-0.0001*y(i)^2-0.0001*x(i)+0.1564*y(i)+0.1325*x(i)*y(i);
    end
    % 更新个体最优解和全局最优解
    for i=1:N
        if z(i) < -64.1290*pbestx(i)^2-0.0001*pbesty(i)^2-0.0001*pbestx(i)+0.1564*pbesty(i)+0.1325*pbestx(i)*pbesty(i)
            pbestx(i) = x(i);
            pbesty(i) = y(i);
        end
        if z(i) < bestval
            bestval = z(i);
            gbestx = x(i);
            gbesty = y(i);
        end
    end
    % 输出迭代结果
    fprintf('第%d代，最优解：%f，坐标（%f，%f）\n',gen,bestval,gbestx,gbesty);
end

增加扰动的惯性加权粒子群算法：

% 增加扰动的惯性加权PSO算法
clear;clc;
N = 50;  % 粒子群数目
maxgen = 100;   % 迭代次数
wmax = 0.9; % 最大惯性权重
wmin = 0.4; % 最小惯性权重
c1 = 2;     % 自我认知学习因子
c2 = 2;     % 社会认知学习因子
r1 = rand(N,1);
r2 = rand(N,1);
Vmax = 10;  % 最大速度
Vmin = -10; % 最小速度
% 初始化种群位置和速度
for i=1:N
    x(i) = rand(1)*100;
    y(i) = rand(1)*100;
    Vx(i) = (rand(1)-0.5)*2*10;
    Vy(i) = (rand(1)-0.5)*2*10;
end
% 计算适应度值
for i=1:N
    z(i) = -64.1290*x(i)^2-0.0001*y(i)^2-0.0001*x(i)+0.1564*y(i)+0.1325*x(i)*y(i);
end
% 初始化全局最优解
[bestval,idx] = min(z);
pbestx = x;
pbesty = y;
gbestx = x(idx);
gbesty = y(idx);
% 迭代寻优
for gen=1:maxgen
    % 计算当前惯性权重
    w = wmax - (wmax-wmin)*gen/maxgen;
    % 更新速度和位置
    for i=1:N
        Vx(i) = w*Vx(i)+c1*r1(i)*(pbestx(i)-x(i))+c2*r2(i)*(gbestx-x(i));
        Vy(i) = w*Vy(i)+c1*r1(i)*(pbesty(i)-y(i))+c2*r2(i)*(gbesty-y(i));
        % 增加扰动
        Vx(i) = Vx(i) + 0.1*randn();
        Vy(i) = Vy(i) + 0.1*randn();
        % 限制速度范围
        Vx(i) = min(max(Vx(i),Vmin),Vmax);
        Vy(i) = min(max(Vy(i),Vmin),Vmax);
        % 更新位置
        x(i) = x(i)+Vx(i);
        y(i) = y(i)+Vy(i);
        % 限制位置范围
        x(i) = min(max(x(i),0),100);
        y(i) = min(max(y(i),0),100);
    end
    % 计算适应度值
    for i=1:N
        z(i) = -64.1290*x(i)^2-0.0001*y(i)^2-0.0001*x(i)+0.1564*y(i)+0.1325*x(i)*y(i);
    end
    % 更新个体最优解和全局最优解
    for i=1:N
        if z(i) < -64.1290*pbestx(i)^2-0.0001*pbesty(i)^2-0.0001*pbestx(i)+0.1564*pbesty(i)+0.1325*pbestx(i)*pbesty(i)
            pbestx(i) = x(i);
            pbesty(i) = y(i);
        end
        if z(i) < bestval
            bestval = z(i);
            gbestx = x(i);
            gbesty = y(i);
        end
    end
    % 输出迭代结果
    fprintf('第%d代，最优解：%f，坐标（%f，%f）\n',gen,bestval,gbestx,gbesty);
end

绘制适应度对比曲线：

% 绘制适应度对比曲线
clear;clc;
N = 50;  % 粒子群数目
maxgen = 100;   % 迭代次数
wmax = 0.9; % 最大惯性权重
wmin = 0.4; % 最小惯性权重
c1 = 2;     % 自我认知学习因子
c2 = 2;     % 社会认知学习因子
r1 = rand(N,1);
r2 = rand(N,1);
Vmax = 10;  % 最大速度
Vmin = -10; % 最小速度
% 初始化种群位置和速度
for i=1:N
    x(i) = rand(1)*100;
    y(i) = rand(1)*100;
    Vx(i) = (rand(1)-0.5)*2*10;
    Vy(i) = (rand(1)-0.5)*2*10;
end
% 计算适应度值
for i=1:N
    z(i) = -64.1290*x(i)^2-0.0001*y(i)^2-0.0001*x(i)+0.1564*y(i)+0.1325*x(i)*y(i);
end
% 初始化全局最优解
[bestval,idx] = min(z);
pbestx = x;
pbesty = y;
gbestx = x(idx);
gbesty = y(idx);
% 传统PSO算法
for gen=1:maxgen
    % 更新速度和位置
    for i=1:N
        Vx(i) = w*Vx(i)+c1*r1(i)*(pbestx(i)-x(i))+c2*r2(i)*(gbestx-x(i));
        Vy(i) = w*Vy(i)+c1*r1(i)*(pbesty(i)-y(i))+c2*r2(i)*(gbesty-y(i));
        % 限制速度范围
        Vx(i) = min(max(Vx(i),Vmin),Vmax);
        Vy(i) = min(max(Vy(i),Vmin),Vmax);
        % 更新位置
        x(i) = x(i)+Vx(i);
        y(i) = y(i)+Vy(i);
        % 限制位置范围
        x(i) = min(max(x(i),0),100);
        y(i) = min(max(y(i),0),100);
    end
    % 计算适应度值
    for i=1:N
        z(i) = -64.1290*x(i)^2-0.0001*y(i)^2-0.0001*x(i)+0.1564*y(i)+0.1325*x(i)*y(i);
    end
    % 更新个体最优解和全局最优解
    for i=1:N
        if z(i) < -64.1290*pbestx(i)^2-0.0001*pbesty(i)^2-0.0001*pbestx(i)+0.1564*pbesty(i)+0.1325*pbestx(i)*pbesty(i)
            pbestx(i) = x(i);
            pbesty(i) = y(i);
        end
        if z(i) < bestval
            bestval = z(i);
            gbestx = x(i);
            gbesty = y(i);
        end
    end
    % 记录适应度值
    f1(gen) = bestval;
end
% 惯性加权PSO算法
for gen=1:maxgen
    % 计算当前惯性权重
    w = wmax - (wmax-wmin)*gen/maxgen;
    % 更新速度和位置
    for i=1:N
        Vx(i) = w*Vx(i)+c1*r1(i)*(pbestx(i)-x(i))+c2*r2(i)*(gbestx-x(i));
        Vy(i) = w*Vy(i)+c1*r1(i)*(pbesty(i)-y(i))+c2*r2(i)*(gbesty-y(i));
        % 限制速度范围
        Vx(i) = min(max(Vx(i),Vmin),Vmax);
        Vy(i) = min(max(Vy(i),Vmin),Vmax);
        % 更新位置
        x(i) = x(i)+Vx(i);
        y(i) = y(i)+Vy(i);
        % 限制位置范围
        x(i) = min(max(x(i),0),100);
        y(i) = min(max(y(i),0),100);
    end
    % 计算适应度值
    for i=1:N
        z(i) = -64.1290*x(i)^2-0.0001*y(i)^2-0.0001*x(i)+0.1564*y(i)+0.1325*x(i)*y(i);
    end
    % 更新个体最优解和全局最优解
    for i=1:N
        if z(i) < -64.1290*pbestx(i)^2-0.0001*pbesty(i)^2-0.0001*pbestx(i)+0.1564*pbesty(i)+0.1325*pbestx(i)*pbesty(i)
            pbestx(i) = x(i);
            pbesty(i) = y(i);
        end
        if z(i) < bestval
            bestval = z(i);
            gbestx = x(i);
            gbesty = y(i);
        end

将以下代码转换为python：function newpop=zmutate(pop,popsize,pm1,pm2,fitness1,M,N,Tn0,Tn1,Q,ST0,maxT,t,maxgen,LCR,ECR,MCR,FC,ICR) %M为辅助坑道数量；N为单元数 x=pop(:,1:2*M+1);%分段点位置 y=pop(:,2*M+2:4*M+2);%是否选择该分段点 z=pop(:,4*M+3:6*M+4);%开挖方向 W=pop(:,6*M+5:8*M+6);%作业班次 lenx=length(x(1,:)); leny=length(y(1,:)); lenz=length(z(1,:)); lenW=length(W(1,:)); avefit=sum(fitness1)/popsize; worstfit=min(fitness1); % sumy=sum(y); % lenz=sumy+1; % lenW=sumy+1; for i=1:popsize %选择popsize次，每次选择一个，输出一个 %随机选择一个染色体 pick=rand; while pick==0 pick=rand; end index=ceil(pick*popsize); f1=fitness1(index); if f1<=avefit % pm=(exp(-t/maxgen))*(pm1-(pm1-pm2)*(f1-avefit)/max(fitness1)-avefit); pm=1/(1+exp(t/maxgen))*(pm1-(pm1-pm2)*(f1-avefit)/max(fitness1)-avefit); else % pm=(exp(-t/maxgen))*pm1; pm=1/(1+exp(t/maxgen))*pm1; end pick=rand; while pick==0 pick=rand; end if pick>pm continue; end % flag0=0; % while(flag0==0) %随机选择变异位置 pick1=rand; pick2=rand; pick3=rand; pick4=rand; while pick1*pick2*pick3*pick4==0 pick1=rand; pick2=rand; pick3=rand; pick4=rand; end posx=ceil(pick1*lenx); posy=ceil(pick2*leny); %x,y变异 randx=randi([1,N-1]); while ismember(randx,x(index,:)) randx=randi([1,N-1]); end b=x(index,posx); x(index,posx)=randx; a=[0 1]; c=y(index,posy); y(index,posy)=setxor(y(index,posy),a); %z,W变异 posz=ceil(pick3*lenz); posW=ceil(pick4*lenW); d=z(index,posz); z(index,posz)=setxor(z(index,posz),a); randW=randi([1,3]); while randW==W(index,posW) randW=randi([1,3]); end e=W(index,posW); W(index,posW)=randW; mpop=[x(index,:),y(index,:),z(index,:),W(index,:)]; mtime=ztime(mpop,M,N,Tn0,Tn1,Q,ST0); mutfit=zcost(mpop,M,N,mtime(:,1),mtime(:,2:2*M+3),mtime(:,2*M+4:2*M+2+N),LCR,ECR,MCR,FC,ICR,Q); if mtime(:,1)>maxT||mutfit<=worstfit x(index,posx)=b; y(index,posy)=c; z(index,posz)=d; W(index,posW)=e; end end newpop=[x,y,z,W]; end

以下是Python代码：

import numpy as np

def zmutate(pop, popsize, pm1, pm2, fitness1, M, N, Tn0, Tn1, Q, ST0, maxT, t, maxgen, LCR, ECR, MCR, FC, ICR):
    # M为辅助坑道数量；N为单元数
    x = pop[:, :2*M+1] # 分段点位置
    y = pop[:, 2*M+2:4*M+2] # 是否选择该分段点
    z = pop[:, 4*M+3:6*M+4] # 开挖方向
    W = pop[:, 6*M+5:8*M+6] # 作业班次

    lenx = x.shape[1]
    leny = y.shape[1]
    lenz = z.shape[1]
    lenW = W.shape[1]

    avefit = np.sum(fitness1) / popsize
    worstfit = np.min(fitness1)

    newpop = np.zeros(pop.shape)

    for i in range(popsize):
        # 随机选择一个染色体
        pick = np.random.rand()
        while pick == 0:
            pick = np.random.rand()
        index = int(np.ceil(pick * popsize))
        f1 = fitness1[index]

        if f1 <= avefit:
            pm = (np.exp(-t / maxgen)) * (pm1 - (pm1 - pm2) * (f1 - avefit) / np.max(fitness1) - avefit)
            pm = 1 / (1 + np.exp(t / maxgen)) * (pm1 - (pm1 - pm2) * (f1 - avefit) / np.max(fitness1) - avefit)
        else:
            pm = (np.exp(-t / maxgen)) * pm1
            pm = 1 / (1 + np.exp(t / maxgen)) * pm1

        pick = np.random.rand()
        while pick == 0:
            pick = np.random.rand()

        if pick > pm:
            newpop[i, :] = pop[index, :]
            continue

        # 随机选择变异位置
        pick1 = np.random.rand()
        pick2 = np.random.rand()
        pick3 = np.random.rand()
        pick4 = np.random.rand()

        while pick1 * pick2 * pick3 * pick4 == 0:
            pick1 = np.random.rand()
            pick2 = np.random.rand()
            pick3 = np.random.rand()
            pick4 = np.random.rand()

        posx = int(np.ceil(pick1 * lenx))
        posy = int(np.ceil(pick2 * leny))

        # x,y变异
        randx = np.random.randint(1, N)
        while randx in x[index, :]:
            randx = np.random.randint(1, N)

        b = x[index, posx]
        x[index, posx] = randx

        c = y[index, posy]
        y[index, posy] = np.setxor1d(y[index, posy], [0, 1])

        # z,W变异
        posz = int(np.ceil(pick3 * lenz))
        posW = int(np.ceil(pick4 * lenW))

        d = z[index, posz]
        z[index, posz] = np.setxor1d(z[index, posz], [0, 1])

        randW = np.random.randint(1, 4)
        while randW == W[index, posW]:
            randW = np.random.randint(1, 4)

        e = W[index, posW]
        W[index, posW] = randW

        mpop = np.concatenate((x[index, :], y[index, :], z[index, :], W[index, :]))
        mtime = ztime(mpop, M, N, Tn0, Tn1, Q, ST0)
        mutfit = zcost(mpop, M, N, mtime[:, 0], mtime[:, 1:2*M+3], mtime[:, 2*M+4:2*M+2+N], LCR, ECR, MCR, FC, ICR, Q)

        if mtime[:, 0] > maxT or mutfit <= worstfit:
            x[index, posx] = b
            y[index, posy] = c
            z[index, posz] = d
            W[index, posW] = e

        newpop[i, :] = np.concatenate((x[index, :], y[index, :], z[index, :], W[index, :]))

    return newpop

注意，Python的数组下标是从0开始的，而MATLAB是从1开始的。此外，Python中的setxor函数的用法与MATLAB有些不同，需要将第二个参数改为一个数组。