Matlab中采用PSO求解整数规划0-1问题

最近做优化方面的算例，发现cplex很难解决整型变量很多的整数规划问题，01变量少的时候还是可以轻松求解，随着01变量数增加求解时间开始指数级增长，在01变量达到100+时直接无法求解了（一直在求解中）...

为了解决这个问题，还是不得不求助智能算法，智能算法的优势在于，虽然给出的结果可能达到局部最优，但肯定会给一个结果。

确定这个方向后，在csdn及github上都搜索了PSO的相关代码，发现主流的有采用01变量及罚函数两种方法，但或者缺少详细原理说明或者采用的并非matlab语言（对于电气学子matlab真的很重要！！）。迫不得已参考了 博主南音小榭PSO算法求解整数规划问题举例的这篇博文，参考里面的代码，将其原本采用的C语言书写的求解代码改为了matlab语言，再次感谢！！！也在这里发出来给大家参考一下，如有不对欢迎指正~

求解的问题：

%% 整数规划0-1
% 目标函数 
    max z = 6x1 + 2x2 + 3x3 + 5x4
% 约束
    4x1 + 2x2 + x3 + 3x4 <= 10
	3x1 - 5x2 + x3 + 6x4 >= 4
	2x1 + x2 + x3 - x4 <= 3
	x1 + 2x2 + 4x3 + 5x4 <= 10
	x[j] = 0 or 1
% 最优解
    z=14,x=(1,0,1,1)

 代码实现：

clc;clear;
%% 初始化
%% PSO参数
person_number = 20;    % 种群个体数目
dimension = 4;         % 搜索变量维度
LDW_ini = 0.9;         % 线性递减权值策略
LDW_end = 0.4;
c1 = 2;c2 = 2;         % 学习因子
maxgen = 100;           % 迭代次数

X_pos = zeros(person_number,dimension);            % 粒子位置变量
V_speed =  zeros(person_number,dimension);         % 粒子速度变量
pbest = zeros(person_number,1);                    % 粒子适应值
max_pbest = 0;                                     % 临时个体最优值
gbest = 0;                                         % 粒子最优值
gbest_temp = 0;                                    % 临时全局最优值——更新历史值时使用
gbest_temp_2 = 0;                                  % 临时全局最优值——约束判断中使用
gbest_index = 1;                                   % 迭代过程中全局最优时粒子位置
gbest_X = zeros(dimension,1);                      % 记录相应的X取值

%% 数据记录
iter_arr = zeros(maxgen,1);                        % 迭代次数
gbest_arr = zeros(maxgen,1);                       % 迭代过程中的最优解

%% 目标约束
const_num = 4;         % 约束数量
C = [6;2;3;5];         % 目标函数系数
A = [4,2,1,3;3,-5,1,6;2,1,1,-1;1,2,4,5];   % 约束系数
B = [10;4;3;10];       % 约束右侧值

%% 主函数
%initial();
[max_pbest,gbest_index,gbest,X_pos,V_speed,pbest] = initial(person_number,dimension,V_speed,pbest,C,X_pos,gbest_index,const_num,A,B);
%gbest = pso(maxgen,LDW_ini,LDW_end,person_number,dimension,max_pbest,gbest_index,gbest,X_pos,V_speed,c1,c2,pbest,C,const_num,A,B,gbest_X,iter_arr,gbest_arr);

iter = 1;
while(iter < maxgen)
   weight = Weight(iter,LDW_ini,LDW_end,maxgen);%Weight(iter);
   % 更新速度和位置
   for i = 1:person_number
      for j = 1:dimension
         r_rand = zeros(2,1);
         for k = 1:2
            r_rand(k,1) = randi([0,9]) / 10;
         end
         V_speed(i,j) = weight * V_speed(i,j) + c1*r_rand(k,1)*(pbest(i,1)-X_pos(i,j)) + c2*r_rand(k,1)*(gbest-X_pos(i,j));
        if V_speed(i,j) < 0 || V_speed(i,j) > 1  % 越界检测，越界后将其限定在0-0.9间
            V_speed(i,j) = randi([0,9]) / 10;
        end
        X_pos(i,j) = X_pos(i,j) + V_speed(i,j);
        if X_pos(i,j) ~= 0 || X_pos(i,j) ~= 1
            X_pos(i,j) = randi([0,1]);
        end
      end
   end
        % 更新个体最优值
        for i = 1:person_number
           pbest(i,1) = Fitness(i,dimension,C,X_pos); 
        end

        max_pbest = 0;
        for i = 1:person_number
           if pbest(i,1) > max_pbest
              max_pbest = pbest(i,1); 
              gbest_index = i;
           else
               max_pbest = max_pbest;
               gbest_index = gbest_index;
           end
        end
        % 约束条件判断
        gbest_temp = Constraint(dimension,gbest_index,const_num,A,B,X_pos,max_pbest);
        % 更新历史最优值
        if gbest_temp > gbest
            gbest = gbest_temp;
            for j = 1:dimension
               gbest_X(j,1) = X_pos(gbest_index,j); 
            end
        else
            gbest = gbest;
            for j = 1:dimension
                gbest_X(j,1) = gbest_X(j,1);
            end
        end
        iter_arr(iter,1) = iter;
        gbest_arr(iter) = gbest;
        iter = iter + 1;
end


%% 计算对应种群的适应度函数
function fit = Fitness(pers, dim,C,X_pos)
    fit = 0;
    for j = 1:dim
        fit = fit + C(j,1) * X_pos(pers,j);
    end
end

%% 判断输入的对应种群数据是否满足约束条件
function gbest_temp_2 = Constraint(dim,index,const_num,A,B,X_pos,max_pbest)
    constraint = zeros(const_num,1);
    gbest_temp_2 = 0;
    for i = 1:const_num
        for j = 1:dim
            constraint(i,1) = constraint(i,1) + A(i,j) * X_pos(index,j);
        end
    end
    if(constraint(1,1) <= B(1,1) && constraint(2,1) >= B(2,1) && constraint(3,1) <= B(3,1) && constraint(4,1) <= B(4,1))
       if max_pbest > gbest_temp_2
           gbest_temp_2 = max_pbest;
       else
           gbest_temp_2 = gbest_temp_2;
       end
    end
end

%% 更新pso权重
function pso_weight = Weight(ITER,LDW_ini,LDW_end,maxgen)
    pso_weight = 0;
    pso_weight = ((LDW_ini - LDW_end) * (maxgen - ITER) / maxgen) + LDW_end;
end

%% 初始化数据
function [max_pbest,gbest_index,gbest,X_pso,V_speed,pbest] = initial(person_number,dimension,V_speed,pbest,C,X_pos,gbest_index,const_num,A,B)
    for i = 1:person_number
       for j = 1:dimension
            X_pso(i,j) = randi([0,1]);             % 初始数据为01变量
            V_speed(i,j) = randi([0,9]) / 10;      % 速度为0-0.9
            % V_speed(i,j) = randi([0,1]);
       end
    end
    for i = 1:person_number
       pbest(i,1) = Fitness(i,dimension,C,X_pos);%(i,dimension);          % 计算每个种群的适应度
    end
    max_pbest = 0;
    for i = 1:person_number
       if pbest(i,1) > max_pbest
          max_pbest = pbest(i,1);                  % 记录初始化种群的适应度最优值
          gbest_index = i;                         % 记录适应度最优时的种群序号
       else
           max_pbest = max_pbest;
           gbest_index = gbest_index;
       end
    end
    gbest = Constraint(dimension,gbest_index,const_num,A,B,X_pos,max_pbest);     % 判断最优种群是否满足约束条件
end