粒子群演算法(PSO)求解订单接受

此代码是我给我朋友写的，如果有人做平行机台订单接受方向或者论文，可以参考此代码，希望对读者有帮助和启示！

可复制直接运行，无需修改！

 

 

/*
PSO algrithm for accept order and schedule
Written by tang jun on February 20, 2020
*/


#include <iostream>
# include <fstream>
#include <iomanip>
#include <math.h>
#include <vector>
#include<random>
#include<ctime>
#include <stdlib.h>


using namespace std;

 


/*
加工时间，收益，延迟惩罚 软时间窗（截止日期前后延长一点），库存成本，序列相依装设，非等效机台
符号索引:
N: 表示订单总数
M: 表示机台总数
i:订单索引，i=1,…,N
m:机台索引，m=1,…,M
参数:
s_0i:订单i在机台上最先被加工所产生的整备时间，i=1,…,N [N,]
s_ij:订单i处理完接着加工订单j在机台上所产生的整备时间，i=1,…,N,j=1,…,N,i≠j。 [N, N]
p_im:订单i在机台m上的加工时间，i=1,…,N。 [N, M]
r_i:订单i的利润，i=1,…,N。 [N,]
w_i:订单i单位延迟处罚权重，i=1,…,N。 [N,]
h_i:订单i单位提早库存成本，i=1,…,N。 [N,]
d_i:订单i的交货时间窗，i=1,…,N。 [ d_i1，d_i2 ] [N , 2]
变数:
C_i:订单i的完工时间，i=1,…,N。 [N,]
T_i:订单i的延迟时间，i=1,…,N。 Ti1= max(0,di1- Ci) >0表示订单i的在时间窗提前完工，提早时间，
有库存成本，否则不考虑 Ti2=max(0,Ci-di2) >0表示订单i 的延迟时间,否则不考虑 [Ti1,Ti2] [N,2]

目标函数：max=ri-Ti1*hi-Ti2*wi i表示被接受的订单

 

*/

 

// 下面为全局变量，便于所有函数或类等都可以使用！

const int N = 20; // total number of orders
const int M = 6; // total number of machines
const int K = 10; // toral number of particle for pso algrithm


// 构建矩阵变量来存储需要的数据，共有2部分：
//第一部分：构建存储随机产生的已知变量

float s_0i[N]; //序列相依装设的时间，用于每台机器第一订单位置
float s_ij[N][N]; //序列相依装设的时间
float p_im[N][M]; // 每台机器对应订单的加工时间
float r_i[N]; //每个订单的利润
float w_i[N]; // 订单i单位延迟处罚权重
float h_i[N]; // 订单i单位提早库存成本
float d_i[N][2]; // 订单i的交货时间窗 [d_i1，d_i2][N, 2]
float C_i[N]; // 订单i的完工时间 [N, ]
float T_i[N][2]; // 订单i的提早和延迟时间


float Randnumber = 10; //随机种子

// PSO算法参数设置

int pso_num_iteration_end = 100000;

float pso_time_end = 100;
float pso_w = 0.2;
float pso_c1 = 0.8;
float pso_c2 = 0.8;
float pso_order_vmin = -8;
float pso_order_vmax = 8;
float pso_machine_vmin = -8;
float pso_machine_vmax = 8;

 

 

// 第二部分：构建所有的粒子订单矩阵及机器矩阵来保存变量

float order_p_gbest[N]; //保存所有粒子中最佳位置
float order_v_gbest[N]; //保存所有粒子中最佳速度
float machine_p_gbest[N]; //保存所有粒子中最佳位置对应的机器最佳位置
float machine_v_gbest[N]; //保存所有粒子中最佳速度对应的机器最佳速度


class util {
// 辅助函数
public:

// 该函数是随机在[LB, UB]中产生浮点数
float RandomRge(float LB, float UB) { float x; x = LB + (UB - LB) / (RAND_MAX) * (float)rand(); return (x); }

// 该函数是随机在[0, 1]中产生浮点数
float RandomF() { float x; x = (float)1.0 / (float)(RAND_MAX) * (float)rand(); return (x); }


// 该函数给定一维数组，输出从大到小的排序的索引值，此函数是一个指针函数，需要建立指针来读取输出数据。
int *sort_seq( float *p, const int size=N) {
float *f;
static int index_int[N];
float temp[N] ;
float temp1 = 0.0;
int index = 0;
int temp_index = 0;
f = p; // 将指针指向了p地址
for (int i = 0; i < size; i++) {
// 该循环将输入的一维变量传值给了temp的变量
temp[i] = *f;
// cout << "temp[i]" << temp[i]<<endl;
f++;
}
for (int i = 0; i < size; i++) {
temp1 = temp[i];
temp_index = i;
for (int j = 0; j < size; j++) {
if ( temp1 > temp[j] ) {
temp1 = temp[j];
temp_index = j;
index = temp_index;
}
else
{
index = temp_index;
}
}
temp[index] = 999999; // 找到对应的索引将其值变一个很大的数
index_int[i] = index; // 将该最小值索引添加到该index_int数组中,最终是输入序列从小到大的索引值
}
return index_int;
}

// 该函数输入类型clock_t的开始时间，就可以计算出从开始时间截止目前时间的耗费时间长,返回是秒
float clock_time(clock_t time_start) {
float time_period;
clock_t time_end;
time_end = clock();
time_period = (time_start - time_end) / CLOCKS_PER_SEC;

return time_period;

}

}U;


void data_product_init() {
float d_i_temp;
srand(Randnumber); // 产生伪随机，在Randnumber不变情况下，以下变量不变
for (int i = 0; i < N; i++) {

s_0i[i] = U.RandomRge(2, 8);

r_i[i] = U.RandomRge(100, 200); //每个订单的利润
w_i[i] = U.RandomRge(2, 10); // 订单i单位延迟处罚权重
h_i[i] = U.RandomRge(2, 10); // 订单i单位提早库存成本

d_i[i][0] = U.RandomRge(15, 50); // 订单i的交货时间窗 [d_i1，d_i2][N, 2]
d_i_temp = U.RandomRge(0, 10); //中间变量，用来加这一部分时间便是截止窗末端
d_i[i][1] = d_i_temp+ d_i[i][0]; // 订单i的交货时间窗 [d_i1，d_i2][N, 2]


C_i[i] = 0; // 根据计算得到，所以不用产生，初始化为0
T_i[i][0] = 0; // 根据计算得到，所以不用产生，初始化为0
T_i[i][1] = 0; // 根据计算得到，所以不用产生，初始化为0

// 对保存全局最佳变量初始化
order_p_gbest[i]=0; //根据计算得到，所以不用产生，初始化为0
order_v_gbest[i]=0; //根据计算得到，所以不用产生，初始化为0
machine_p_gbest[i]=0; //根据计算得到，所以不用产生，初始化为0
machine_v_gbest[i]=0; //根据计算得到，所以不用产生，初始化为0


for (int j = 0; j < N; j++) {
s_ij[i][j]= U.RandomRge(2, 8); //序列相依装设的时间
}


for (int m = 0; m < M; m++) {
p_im[i][m]= U.RandomRge(5, 20); // 每台机器对应订单的加工时间
}

}


}


class PSO {

public:
//创建每个粒子的变量
float order_position_seq[N]; //创建并保存订单位置序列的变量 16
float order_speed_seq[N]; //创建并保存订单速度序列的变量
float order_p_best[N]; //保存该粒子订单位置最佳位置
float order_v_best[N]; //保存该粒子订单最佳速度
// 创建机器变量，随后迭代保存最佳需要对应粒子保存的最佳
float machine_position_seq[N]; //创建机器位置序列的变量 15
float machine_speed_seq[N]; //创建机器速度序列的变量
float machine_p_best[N]; //保存该粒子机器最佳位置
float machine_v_best[N]; //保存该粒子机器最佳速度
PSO() {
// 构造函数用来随机初始化变量
for (int i = 0; i < N; i++) {
order_position_seq[i]= U.RandomRge(1,N);
order_speed_seq[i] = U.RandomRge(1, N);
order_p_best[i]= order_position_seq[i];
order_v_best[i]= order_speed_seq[i];
machine_position_seq[i] = U.RandomRge(0, M+1-0.00001);
machine_speed_seq[i] = U.RandomRge(0, M + 1 - 0.00001);
machine_p_best[i]= machine_position_seq[i];
machine_v_best[i]= machine_speed_seq[i];
}
}

}particle[K]; //particle[K] 表示已经对象化了，产生了K个粒子


float object(float *order_seq,float *machine_seq) {
// 构建目标函数
float *order; //建立浮点型指针用来读取order_seq的数据
float *machine;//建立浮点型指针用来读取machine_seq的数据
int *index; //建立订单排序函数指针，用来读取排序函数的输出值
float order_sequence[N]; //用来保存order序列
int order_index[N];
float machine_sequence[N];//用来保存机台序列
order = order_seq; //该指针指向order序列
machine = machine_seq;//该指针指向machine序列
for (int i = 0; i < N; i++) {
// 该循环将值赋给上面的变量。
order_sequence[i] = *order;
machine_sequence[i] = *machine;
order++;
machine++;
}

index = U.sort_seq(order_sequence); //用index指针指向输出值的首地址

int machine_sequence_int[N] = { 0 };

for (int i = 0; i < N; i++) {
// 将订单序列从小到大的索引已经赋值给order_index中了。
order_index[i] = *(index+i);
machine_sequence_int[i] = int(machine_sequence[i]);
}


int machine_count[M]; // 保存每个机台对应加工订单的数量
for (int m = 0; m < M; m++) { machine_count[m] = 0; } // 赋值为0

 

for (int i = 0; i < N; i++ ) {
for (int m = 0; m < M; m++) {
if (machine_sequence_int[i]==m+1) {
// 利用m+1可以排除机台为0的不用计算进去
machine_count[m]++;
}
}
}
// 下面较为复杂，按照每台机器取计算目标值
int order_temp1 = 0;
int order_temp_i = -1;
int order_temp_j = -1;
float object_result = 0.0;

int temp_cont = 0; // 使用计数方法来处理机台的第一订单

for (int m = 0; m < M; m++) {
if (machine_count[m] == 1) {
// 该条件说明机台只有一个订单加工，则只有装设准备时间需要s_0i
for (int i = 0; i < N; i++) {
order_temp1 = order_index[i]; // 将第i个位置加工的订单索引给变量order_temp1
if (machine_sequence_int[order_temp1] == m + 1) {
// 对应机台是否等于该循环的机台，如果是就记住该索引号，停止循环，因为只有一个订单
order_temp_i = order_temp1;
break; // 跳出当前循环
}
}

C_i[order_temp_i]=s_0i[order_temp_i]+ p_im[order_temp_i][machine_sequence_int[order_temp_i]];
T_i[order_temp_i][0]=d_i[order_temp_i][0]- C_i[order_temp_i]; // 如果大于0表示提前完成有库存成本
T_i[order_temp_i][1] = C_i[order_temp_i]-d_i[order_temp_i][1]; // 如果大于0表示延迟完成有惩罚权重
if (T_i[order_temp_i][0] < 0) { T_i[order_temp_i][0] = 0; }
if (T_i[order_temp_i][1] < 0) { T_i[order_temp_i][1] = 0; }
// 计算目标值 ri-Ti1*hi-Ti2*wi
object_result = object_result + r_i[order_temp_i] - h_i[order_temp_i] * T_i[order_temp_i][0] - w_i[order_temp_i] * T_i[order_temp_i][1];

}

if (machine_count[m] > 1) {
// 该条件说明机第（m+1)的机台有多个订单加工
temp_cont = 0; // 将其赋值为0
for (int i = 0; i < N; i++) {
order_temp1 = order_index[i]; // 将第i个位置加工的订单索引给变量order_temp1
if (machine_sequence_int[order_temp1] == m + 1) {
// 对应机台是否等于该循环的机台，如果是就记住该索引号
temp_cont++;//计数，若为1表示该机台首次加工该订单，需要s0i
if (temp_cont==1) {
order_temp_i = order_temp1;
C_i[order_temp_i] = s_0i[order_temp_i] + p_im[order_temp_i][machine_sequence_int[order_temp_i]];
T_i[order_temp_i][0] = d_i[order_temp_i][0] - C_i[order_temp_i]; // 如果大于0表示提前完成有库存成本
T_i[order_temp_i][1] = C_i[order_temp_i] - d_i[order_temp_i][1]; // 如果大于0表示延迟完成有惩罚权重
if (T_i[order_temp_i][0] < 0) { T_i[order_temp_i][0] = 0; }
if (T_i[order_temp_i][1] < 0) { T_i[order_temp_i][1] = 0; }
// 计算目标值 ri-Ti1*hi-Ti2*wi
object_result = object_result + r_i[order_temp_i] - h_i[order_temp_i] * T_i[order_temp_i][0] - w_i[order_temp_i] * T_i[order_temp_i][1];

}
else {
order_temp_j = order_temp1;
C_i[order_temp_j] = C_i[order_temp_i]+s_ij[order_temp_i][order_temp_j] + p_im[order_temp_j][machine_sequence_int[order_temp_j]];
T_i[order_temp_j][0] = d_i[order_temp_j][0] - C_i[order_temp_j]; // 如果大于0表示提前完成有库存成本
T_i[order_temp_j][1] = C_i[order_temp_j] - d_i[order_temp_j][1]; // 如果大于0表示延迟完成有惩罚权重
if (T_i[order_temp_j][0] < 0) { T_i[order_temp_j][0] = 0; }
if (T_i[order_temp_j][1] < 0) { T_i[order_temp_j][1] = 0; }
// 计算目标值 ri-Ti1*hi-Ti2*wi
object_result = object_result + r_i[order_temp_j] - h_i[order_temp_j] * T_i[order_temp_j][0] - w_i[order_temp_j] * T_i[order_temp_j][1];
order_temp_i = order_temp_j; // 将处理订单传给上一个订单变量，下次变成了这个订单
}

}
}

}

}
return object_result;

}


float pso_run() {

// 对生产变量初始化,直接调用函数
data_product_init();

 

 


//下面循环将是在k个粒子中找到全局最佳解将其值传给
float object_old=0.0;
float object_new=0.0;
float object_invidual = 0.0;

//在K个粒子中寻找最佳解，并将其赋值给最佳值变量中
for (int k = 0; k < K; k++) {

if (k == 0) {
object_old = object(particle[k].order_p_best, particle[k].machine_p_best);
for (int i = 0; i < N; i++) {
order_p_gbest[i] = particle[k].order_p_best[i];
order_v_gbest[i] = particle[k].order_v_best[i];
machine_p_gbest[i] = particle[k].machine_p_best[i];
machine_v_gbest[i] = particle[k].machine_v_best[i];
}
}
else {
object_new = object(particle[k].order_p_best, particle[k].machine_p_best);
if (object_new > object_old) {
object_old = object_new; // 如果大就换值
for (int i = 0; i < N; i++) {
order_p_gbest[i] = particle[k].order_p_best[i];
order_v_gbest[i] = particle[k].order_v_best[i];
machine_p_gbest[i] = particle[k].machine_p_best[i];
machine_v_gbest[i] = particle[k].machine_v_best[i];
}
}
}

//cout << "object=" << object_old << endl;

}

 


float order_p_new[N] = { 0 };
float order_v_new[N] = { 0 };
float machine_p_new[N] = { 0 };
float machine_v_new[N] = { 0 };

int iteration_count = 0;
clock_t time_start;
float iteration_time = 0;

time_start = clock(); //建立停止时间
while (iteration_count < pso_num_iteration_end) {

if (iteration_time > pso_time_end) { break; }

// 下面是对每只粒子进行解的迭代
for (int k = 0; k < K; k++) {
// 对粒子循环
for (int i = 0; i < N; i++) {
// 对订单及机台进行领域解的更新
//对订单的解进行更新
order_v_new[i] = pso_w * particle[k].order_position_seq[i] +
pso_c1 * U.RandomF() * (particle[k].order_p_best[i] - particle[k].order_position_seq[i]) +
pso_c2 * U.RandomF() * (order_p_gbest[i] - particle[k].order_position_seq[i]);
// 对速度进行一个判断
if (order_v_new[i] > pso_order_vmax) { order_v_new[i] = pso_order_vmax; }
if (order_v_new[i] < pso_order_vmin) { order_v_new[i] = pso_order_vmin; }
order_p_new[i] = particle[k].order_position_seq[i] + order_v_new[i];
//对机台的解进行更新
machine_v_new[i] = pso_w * particle[k].machine_position_seq[i] +
pso_c1 * U.RandomF() * (particle[k].machine_p_best[i] - particle[k].machine_position_seq[i]) +
pso_c2 * U.RandomF() * (machine_p_gbest[i] - particle[k].machine_position_seq[i]);
// 对速度进行一个判断
if (machine_v_new[i] > pso_machine_vmax) { machine_v_new[i] = pso_machine_vmax; }
if (machine_v_new[i] < pso_machine_vmin) { machine_v_new[i] = pso_machine_vmin; }
machine_p_new[i] = particle[k].machine_position_seq[i] + machine_v_new[i];
// 对机器做一些限制,因为超出机台数量那就有问题了
if (machine_p_new[i] > M) { machine_p_new[i] = M; }
if (machine_p_new[i] < 0) { machine_p_new[i] = 0; }
}
// 求解未更新的解与已经更新的解
object_old = object(particle[k].order_position_seq, particle[k].machine_position_seq);
object_new = object(order_p_new, machine_p_new);

// 下面的代码是根据object_old与object_new是否需要更换解
if (object_new >= object_old) {
object_old = object(order_p_gbest, machine_p_gbest); // 求解全局最佳解的目标值
// 满足该条件表示需要更换该解，我们先把个体最佳解给更新了
for (int i = 0; i < N; i++) {
particle[k].order_position_seq[i]=order_p_new[i];
particle[k].order_speed_seq[i]=order_v_new[i];
particle[k].machine_position_seq[i]=machine_p_new[i];
particle[k].machine_speed_seq[i]=machine_v_new[i];

object_invidual= object(particle[k].order_p_best, particle[k].machine_p_best); // 求解当前粒子个体最佳解的目标值
if (object_invidual> object_new) {
particle[k].order_p_best[i] = order_p_new[i];
particle[k].order_v_best[i] = order_v_new[i];
particle[k].machine_p_best[i] = machine_p_new[i];
particle[k].machine_v_best[i] = machine_v_new[i];

}

if (object_new > object_old) {
// 该条件满足需要更新全局最佳解
order_p_gbest[i] = particle[k].order_position_seq[i];
order_v_gbest[i] = particle[k].order_speed_seq[i];
machine_p_gbest[i] = particle[k].machine_position_seq[i];
machine_v_gbest[i] = particle[k].machine_speed_seq[i];
}
}
}
}

iteration_count++; //迭代次数叠加
iteration_time = U.clock_time(time_start);

float object_end = 0.0; // 打印出最佳值
object_end = object(order_p_gbest, machine_p_gbest);
cout << "value iteratived by pso algrithm :" << object_end << endl;

}

 

 

float object_end = 0.0; // 打印出最佳值
object_end= object(order_p_gbest, machine_p_gbest);
cout << "value finished pso algrithm :" << object_end<<endl;

return object_end;
}

 

 

 


int main(){


pso_run();


while (0);

return 0;
}