pso.cpp

改进的pso算法：引入变异算子

#include <cmath>
#include <iostream>
#include <fstream>
#include <ctime>
using namespace std;
const bool m_bUseSeed = 1;		//是否使用种子
const int PNum = 16;			//处理器的个数，即n
const int iAgentDim = 128;		//任务的个数，即m
const int iRangL = 0;			//个体的取值范围，最小为0
const int iRangR = PNum-1;		//最大为处理器个数－1
const int iPSONum=20;			//粒子数
int gbest[iAgentDim];			//全局最优值的坐标
/*PSO参数设定 */
int iCurStep=0;				//当前迭代次数
const int iStep=1000;		//最大迭代次数
//--下面的值，要具体程序中具体的修改，根据你优化的函数来修改--//
const double w_s=0.9;		//起始惯性系数
const double w_e=0.4;		//中止惯性系数
double w=0.9;				//惯性系数
const double c1=2;
const double c2=2;
int EXT[iAgentDim][PNum];				//记录任务机器分配表的矩阵
const int vmax=abs(iRangR-iRangL)/2;	//最大飞行速度
const int vmin=1;						//最小飞行速度
const int TRANSGENEDEPTH=4;				//转基因深度
int EXCGENE[iAgentDim][PNum];			//存储优良基因的矩阵
const double S0=0.8;					//个体相似度阈值
const double C0=0.8;					//群体多样性阈值		
ifstream fin;
ofstream fout;
#define rnd(low,uper)((double)(rand()/(double)RAND_MAX)*((double)(uper)-(low))+(low))//这个东西，返回low ,uper之间的一个值
void CopyMatrix(int srcM[iAgentDim][PNum], int desM[iAgentDim][PNum], int m, int n)//将m*n的矩阵srcM复制给desM
{
	for (int i=0;i<m;i++)
		for (int j=0;j<n;j++)
		{
			desM[i][j] = srcM[i][j];
		}
}

void CopyVector(int* srcV, int* desV, int n)	//将n维向量srcV复制给desV
{
	for(int i=0;i<n;i++)
		desV[i] = srcV[i];
}

class Agent					//单个的粒子，也就是一只鸟 
{
private:
	int tMachineId;
	int tMaxLength;
public:
	int dpos[iAgentDim];	//位置，也就是各个任务分配到哪个处理器
	int dpbest[iAgentDim];	//维护一个"自己"找到的最优分配的方案
	int dv[iAgentDim];		//速度，限定为整数
	int m_dFitness;			//粒子在当前位置的调度长度
	int m_dBestfitness;		//粒子历史上的最短调度长度
	int m_pMachineId;			//该个体对应的问题机器号
	Agent()					//初始化
	{ 
		srand( (unsigned)(time( NULL )+rand()) );
		int i=0;
		for(i=0;i<iAgentDim;i++)
		{
			dpos[i]=rnd(iRangL,iRangR);	//初始化位置
			dpbest[i]=dpos[i];	//初始化最佳位置
			dv[i] = abs(iRangL-iRangR)/2;			//初始化速度
		}
	}
	//判断和个体t是否相似
	bool AsSame(Agent &t)
	{
		double count = 0;
		for (int i=0;i<iAgentDim;i++)
	// 		if (a.dpos[i]==b.dpos[i])
			if (abs(dpos[i]-t.dpos[i])<1)
				count++;
		if ((count/iAgentDim)>=S0)
			return true;
		else return false;
	}
	int CalculateFitness()
	{
		//求该方案的调度长度
		int length = 0;		//每台机器调度长度
		int maxLength = 0;	//所有机器最大调度长度，即这个方案的调度长度
		for (int i=0;i<PNum;i++)	//依次求每台机器的调度长度
		{
			length = 0;
			for (int j=0;j<iAgentDim;j++)	
			{
				if (i==dpos[j])		//若任务j由机器i处理
					length += EXT[j][dpos[j]];
			}
			if (length>maxLength)
			{
				maxLength = length;
				tMachineId = i;		//问题机器
			}
		}
		tMaxLength = maxLength;
		return maxLength;
	}
	void UpdateFitness()
	/*calculate the fitness and find out the best fitness,record*/
	{
		m_dFitness = CalculateFitness();
		//找到一个更好的值后，更新 m_dBestfitness
		if (m_dFitness<m_dBestfitness)
		{
			m_dBestfitness=m_dFitness;
			int i=0;
			for(i=0;i<iAgentDim;i++)
			{
				dpbest[i]=dpos[i];
			}
		}
		m_pMachineId = tMachineId;	//更新问题机器
	}
	void UpdatePos()		//粒子移动
	{
		for(int i=0;i<iAgentDim;i++)
		{
			if ( w_s&&w_e!=0 )
				w = (w_e-w_s)*iCurStep/iStep + w_s;
			dv[i]=w*dv[i]+c1*rnd(0,1)*(dpbest[i]-dpos[i])+c2*rnd(0,1)*(gbest[i]-dpos[i]);
			//限制飞行速度不能超过vmax，不能小于vmin	
			if (dv[i]>vmax && dv[i]>0)
				dv[i] = vmax;
			else if (dv[i]<-vmax && dv[i]<0)
				dv[i] = -vmax;
			if (dv[i]>-vmin && dv[i]<0)
				dv[i] = -vmin;
			else if (dv[i]<vmin && dv[i]>0)
				dv[i] = vmin;
			dpos[i]+=dv[i];			//X(n+1) = X(n)+V(n)
			if (dpos[i]>PNum-1)
				dpos[i] = PNum-1;
			if (dpos[i]<0)
				dpos[i] = 0;
		}
	}
	/*转基因算子,excellentGene为优良基因,transDepth为转基因深度度*/
	void Transgene(int excellentGene[iAgentDim][PNum], int transDepth)//O(transDepth*iAgentDim*iAgentDim)
	{
		int oldGene[iAgentDim];
		CopyVector(dpos, oldGene, iAgentDim);	//保留原来的基因
		int oldFitness = m_dFitness;			//保留原来的适应值
		int newFitness = 0;						//改进后的适应值
		int i,j,k;
/************************************************************************/
/* 从优良基因中接种                                                     */
/************************************************************************/
		bool flag = false;						//标记是否找到改进
		int pMachine, aidMachine;
		for (i=0;i<transDepth;i++)
		{
			pMachine = m_pMachineId;//找到问题机器
			for (j=0;j<iAgentDim;j++)	//对每一个任务
			{
				if (dpos[j]==pMachine)	//若该任务在问题机器上
				{
					//获得目标机器
					aidMachine = excellentGene[j][transDepth];
// 					if (aidMachine==dpos[j]) break;
					dpos[j] = aidMachine;
					newFitness = CalculateFitness();
					if (newFitness<oldFitness)			//若有改进则退出
					{
						UpdateFitness();
						flag = true;
						break;
					}
					//若无改进，则进行逐位转基因操作
					dpos[j] = oldGene[j];	//恢复原值
					//对每一个任务，若它在目标机器上，则与任务j交换		
					for (k=0;k<iAgentDim;k++)					
					{
						if (dpos[k]==aidMachine)
						{
							int tmp = dpos[j];//交换两任务
							dpos[j] = aidMachine;
							dpos[k] = tmp;
							newFitness = CalculateFitness();
							if (newFitness<oldFitness)			//若有改进则退出
							{
								UpdateFitness();
								flag = true;
								break;
							}
						}
					}
				}
				if (flag==true) break;
			}
			if (flag==true) break;
		}
		if (oldFitness==m_dFitness)	//若没有改进，恢复原来的状态
			CopyVector(oldGene, dpos, iAgentDim);
/************************************************************************/
/* 从优良个体中接种                                                     */
/************************************************************************/		
// 		CopyVector(dpos, oldGene, iAgentDim);	//保留原来的基因
// 		oldFitness = m_dFitness;			//保留原来的适应值
// 		int gBegin = (int)rnd(0,127);
// 		int gEnd = gBegin+iAgentDim/4;
// 		if (gEnd>iAgentDim) gEnd = iAgentDim;
// 		for (i=gBegin;i<gEnd;i++)
// 		{
// 			dpos[i] = gbest[i];
// 		}
// 		newFitness = CalculateFitness();
// 		if (newFitness<oldFitness)			//若有改进则更新
// 			UpdateFitness();
// 		if (oldFitness==m_dFitness)	//若没有改进，恢复原来的状态
// 			CopyVector(oldGene, dpos, iAgentDim);
	}
	void Mutate(double prob)//O(20*transDepth-1*iAgentDim*iAgentDim)
	{
		srand( (unsigned)(time( NULL )+rand()) ); 
		for (int i=0;i<iAgentDim;i++)
		{
			if (prob>=rnd(0,1))		//以概率prob变异
				dpos[i] = (int)rnd(0, PNum);
		}
		for (i=0;i<20;i++)//O(20*transDepth-1*iAgentDim*iAgentDim)
			Transgene(EXCGENE, PNum-1);	//变异以后进行深度转基因操作	
	}

};

class PSO	//这是粒子群，也就是鸟群了
{
private:
	Agent agents[iPSONum];
	int m_iTempPos;
	int m_iBadAgent;	//记录最差个体
	int m_iBadFitness;	//最差适应值
	int caled[iPSONum];	//记录相似个体
public:
	int m_dBestFitness;		//鸟群找到的最优值，即最短调度长度
	void Init();
	void Search();
	bool TheSame(int a, int b);	//判断个体a,b是否相似
	double CalComp();				//计算种群相似性
	int maxAgentId;					//浓度最大的个体
	double maxAgentChro;				//浓度最大个体的浓度
	void Mutate(int id, double prob);//对与个体id相似的所有个体进行概率为prob 的变异
};

bool PSO::TheSame(int a, int b)
{
	double count = 0;
	for (int i=0;i<iAgentDim;i++)
	{
		if (agents[a].dpos[i]==agents[b].dpos[i])
			count++;
	}
	if ((count/iAgentDim)>=S0)
		return true;
	else return false;
}
//计算个体浓度
double PSO::CalComp()
{
	int i,j;
	double count=0;
	double chroma=0;//浓度
	int id=0;
	for (i=0;i<iPSONum;i++)
		caled[i] = -1;
	for (i=0;i<iPSONum;i++)	//依次求和i个体相似的个体数
	{
		count = 0;
		if (caled[i]==-1)
		{
			caled[i] = i;	//个体i本身已计算
			count = 1;
		}
		for (j=i+1;j<iPSONum;j++)
		{
			if ((-1==caled[j]) && TheSame(i, j))
			{
				count++;
				caled[j] = i;//表示已计算过j个体
			}
		}
		if (chroma<count/iPSONum)	//更新浓度
		{	
			id = i;			//记录最浓个体
			chroma = count/iPSONum;
		}
	}
	maxAgentId = id;
	maxAgentChro = chroma;
	return chroma;
}

void PSO::Mutate(int id, double prob)//对与个体id相似的所有个体进行概率为prob 的变异
{
// 	for (int i=0;i<iPSONum;i++)
// 		if (TheSame(agents[id], agents[i]))
// 			agents[i].Mutate(prob);
}

void PSO::Search()
{
	iCurStep = 0;
	while( iCurStep<iStep )//最大迭代次数//迭代最大为O(iPSONum*iAgentDim*PNum,20*128*16)
		//so,O(iStep*iPSONum*iAgentDim*PNum,1000*20*128*16)
	{
		m_iTempPos=INT_MAX;
		int i;
		for(i=0;i<iPSONum;i++)  //找鸟群中有没有更好的解，如果有，记录下来//O(iPSONum)
		{
			if ( agents[i].m_dBestfitness<m_dBestFitness ) 
			{
				m_dBestFitness=agents[i].m_dBestfitness;
				m_iTempPos=i;	//找到到的最好解的位置
			}
		}
		if (m_iTempPos!=INT_MAX)
		{
			int j;
			for(j=0;j<iAgentDim;j++)//任务的个数，即m
			{
				gbest[j]=agents[m_iTempPos].dpos[j];//记录全局最优解的各个坐标
			}
		}
		//printf("The best is %f \n",m_dBestFitness);
		//下一次迭代
		for(i=0;i<iPSONum;i++)//O(iPSONum*iAgentDim*PNum,20*128*16)
		{
			agents[i].UpdatePos();//O(iAgentDim,128)
			agents[i].UpdateFitness();//O(iAgentDim*PNum,128*16)
		}   
		for(i=0;i<iPSONum;i++)  //找群体中最差的解，进行转基因操作
		{
			m_iBadFitness = 0;
			if ( agents[i].m_dBestfitness>m_iBadFitness ) 
			{
				m_iBadFitness = agents[i].m_dBestfitness;
				m_iBadAgent = i;	//找到最差个体
			}				
		}
		agents[m_iBadAgent].Mutate(1);//对最差个体进行变异//O(iAgentDim,128)
		for(i=0;i<iPSONum;i++)  //对所有个体进行转基因操作O(iPSONum*transDepth*iAgentDim*iAgentDim)为转基因深度度
		{
			agents[i].Transgene(EXCGENE, TRANSGENEDEPTH);//O(transDepth*iAgentDim*iAgentDim)为转基因深度度
		}
// 		agents[m_iBadAgent].Transgene(EXCGENE, PNum);//对最差个体进行深度转基因操作
		
		if (iCurStep%100==0)
		{
			printf("%d-", iCurStep);
			printf("%d-", m_dBestFitness);
			printf("%.3f ", CalComp());//O(iPSONum*iPSONum/100)
			if (CalComp()>C0)//O(iPSONum*iPSONum/100*iPSONum*20*PNum-1*iAgentDim*iAgentDim)
			{
				for (i=0;i<iPSONum;i++)
				{
					if (caled[i]==maxAgentId)
						agents[i].Mutate(maxAgentChro);//O(20*transDepth-1*iAgentDim*iAgentDim)
				}
			}
		}
		iCurStep++;
	}
	printf("The best result is: %d after %d steps. \n", m_dBestFitness, iCurStep);
	{
		for (int i=0;i<iAgentDim;i++)   
			printf(" %d ",gbest[i]);
		printf("\n");
	}
}

void PSO::Init()//初始化，
{
	int i=0;
	m_dBestFitness=INT_MAX;
	srand( (unsigned)(time( NULL )+rand()) ); 
	for(i=0;i<iPSONum;i++)//粒子数////O(iPSONum*iAgentDim*PNum)
	{ 
		agents[i].m_dBestfitness =INT_MAX;//将m_dBestfitness赋值为一个大的值，目的是找最小值，
		agents[i].UpdateFitness();//O(iAgentDim*PNum,128*16)
	}
	if (m_bUseSeed)		//若使用种子
	{
		ifstream fin("SMM's result.txt");
		cout<<"种子是： "<<endl;
		for (i=0;i<iAgentDim;i++)//任务的个数，即m
		{
			fin>>agents[0].dpos[i];
			cout<<agents[0].dpos[i]<<" ";
		}
		fin.close();
		cout<<endl;
		agents[0].UpdateFitness();
	}
}

void main()
{
	fin.open("input.txt");
	fout.open("output.txt");
	int i, j, k;
	clock_t tstart, tend;
	tstart = clock();
	for (i=0;i<iAgentDim;i++)//任务个数
		for (j=0;j<PNum;j++)//处理器的个数，即n
		{
			fin>>EXT[i][j];
		}
//	printf("The last num is: %d 。 \n", EXT[iAgentDim-1][PNum-1]);
	//求解优良基因
	int tmpEXT[iAgentDim][PNum];
	CopyMatrix(EXT, tmpEXT, iAgentDim, PNum);//O(iAgentDim*PNum)
	for (i=0;i<iAgentDim;i++)//O(iAgentDim*PNum*PNum)
	{
		int ibest = 0, tlength = INT_MAX;
		for (j=0;j<PNum;j++)
		{
			tlength = INT_MAX;
			for (k=0;k<PNum;k++)
			{
				if (tmpEXT[i][j]<tlength)
				{
					ibest = j;
					tlength = tmpEXT[i][j];
				}
			}
			tmpEXT[i][ibest] = INT_MAX;
			EXCGENE[i][j] = ibest;		//求得第i好的机器号为ibest，保存
		}
	}

	int best = INT_MAX;
	for (i=0;i<1;i++)
	{
		PSO pso;
		pso.Init();
		pso.Search();
		printf("The %d result is: %d 。 \n", i, pso.m_dBestFitness);
		if (best>pso.m_dBestFitness)
			best = pso.m_dBestFitness;
	}
	printf("The average best result is: %d 。 \n", best);
	tend = clock();
	cout<<"总时间:"<<tend-tstart<<endl;
}