pso_learning.cpp

粒子群优化算法

#include "stdafx.h"
#include "stdio.h"
#include "stdlib.h"
#include "time.h"
#include "math.h"
#include "nrutil.h"
#include "hmm.h"
#include "pso_learning.h"



//全局变量定义

double W=1.0;
double C1=1.8;
double C2=1.8;
double VMAX=1;
double XMIN=0.0;
double XMAX=1.0;
double P[DIMENSION]; //全局最优解
double PBEST;        //全局最优值  种群历史最优


struct indi
{
	double number[DIMENSION];
	double best[DIMENSION];
	double bestfitness;//粒子历史最优
	double fitness;
	double speed[DIMENSION];
}individual[POPSIZE];




void pso(HMM *phmm, int **O, int LearnN, int *T)
{
	int i,k,t,total=0;
	float sum=0;
	initiate(O, LearnN, T);
	for(i=0;i<MaxGen;i++)//500代演化
	{
		W=float(1.0-i*0.6/499);
		for(k=0;k<POPSIZE;k++)
		{
			for(t=0;t<DIMENSION;t++)
			{
				individual[k].speed[t]=W*individual[k].speed[t]
							           +C1*rdft()*(individual[k].best[t]-individual[k].number[t])
									   +C2*rdft()*(P[t]-individual[k].number[t]);
				
				if(individual[k].speed[t]>VMAX)
					individual[k].speed[t]=VMAX;

				individual[k].number[t]=individual[k].number[t]+individual[k].speed[t];
                
                if(individual[k].number[t]<XMIN)
					individual[k].number[t]=2*XMIN-individual[k].number[t];
				if(individual[k].speed[t]>XMAX)
					individual[k].number[t]=2*XMAX-individual[k].number[t];
				FV(k);
			}
			calculation(k,O,LearnN,T);
			localbest(k);
		}
		globalbest(1);

	}
	
    //********************************************************
    int s1;

	for(i=0;i<DIMENSION;i++)
	{
		if(i<phmm->N*phmm->N)
		{
			phmm->A[int((i+1)/phmm->N)][(i+1)%phmm->N]=P[i];
		}
		else if(i>=phmm->N*phmm->N && i<phmm->N*(phmm->N + phmm->M))
		{
			s1=i-phmm->N*phmm->N+1;
			phmm->B[int(s1/phmm->N)][s1%phmm->M]=P[i];
		}
		else
		{
			s1=i-phmm->N*(phmm->N+phmm->M)+1;
			phmm->pi[s1]=P[i];
		}
	}
	
}

//程序初始化定义
void initiate(int **O, int LearnN, int *T)
{
	int i,j;
	double sum1=0,sum2=0,sum3=0;

	for(i=0;i<POPSIZE;i++)
	{
		for(j=0;j<DIMENSION;j++)
		{
			individual[i].number[j]=rdft()*(XMAX-XMIN)+XMIN;
			if(j<pow(HmmN,2))
				sum1=sum1+individual[i].number[j];
			else if(j>(pow(HmmN,2)-1)&&j<HmmN*(HmmN+HmmM))
				sum2=sum2+individual[i].number[j];
			else
				sum3=sum3+individual[i].number[j];
		}
		for(j=0;j<DIMENSION;j++)
		{
			if(j<pow(HmmN,2))
				individual[i].number[j]=individual[i].number[j]/sum1;
			else if(j>(pow(HmmN,2)-1)&&j<HmmN*(HmmN+HmmM))
				individual[i].number[j]=individual[i].number[j]/sum2;
			else
				individual[i].number[j]=individual[i].number[j]/sum3;
		}

	}


	for(i=0;i<POPSIZE;i++)
		for(j=0;j<DIMENSION;j++)
			individual[i].speed[j]=VMAX*rdft();
    for(i=0;i<POPSIZE;i++)
		for(j=0;j<DIMENSION;j++)
			individual[i].best[j]=individual[i].number[j];
	for(i=0;i<POPSIZE;i++)
	{
		calculation(i, O, LearnN, T);//计算适应函数值
        individual[i].bestfitness=individual[i].fitness;
	}
	globalbest(0);
}

//粒子历史最优
void localbest(int number)
{
	int i;
	if(individual[number].bestfitness>individual[number].fitness)
	{
		for(i=0;i<DIMENSION;i++)
		{
			individual[number].best[i]=individual[number].number[i];
		}
       individual[number].bestfitness=individual[number].fitness;
	}
}

//种群历史最优
void globalbest(int number)
{
	int i,j;
	double s=0;
	int flag=0;
	if(number==0)
	{
		s=individual[0].fitness;
		flag=0;
		for(i=1;i<POPSIZE;i++)
			if(individual[i].fitness<s)
			{ 
				s=individual[i].fitness;
				flag=i;
			}
         for(i=0;i<DIMENSION;i++)
			 P[i]=individual[flag].number[i];
		 PBEST=individual[flag].fitness;
	}
	else
	{
		for(i=0;i<POPSIZE;i++)
			if(individual[i].bestfitness<PBEST)
			{
				for(j=0;j<DIMENSION;j++)
					P[j]=individual[i].best[j];
				PBEST=individual[i].bestfitness;
			}
	}
}

//计算适应值
void calculation(int num, int **O, int LearnN, int *T)
{
	//**********************************************************
	int i, j, k,t;
	HMM *phmm;
	phmm->M = HmmM;
 	phmm->N = HmmN;
 	phmm->A = (double **) dmatrix(1, phmm->N, 1, phmm->N);
	for (i = 1; i <= phmm->N; i++) 
	{
		for (j = 1; j <= phmm->N; j++) 
		{
			phmm->A[i][j] = individual[num].number[(i-1)*phmm->N+j-1]; 
		}
	}
 
	phmm->B = (double **) dmatrix(1, phmm->N, 1, phmm->M);

	for (j = 1; j <= phmm->N; j++) 
	{
		for (k = 1; k <= phmm->M; k++) 
		{
			phmm->B[j][k] = individual[num].number[phmm->N*phmm->N+(j-1)*phmm->N+k-1]; 
		}
	}
 
	phmm->pi = (double *) dvector(1, phmm->N);
	for (i = 1; i <= phmm->N; i++) 
	{
		phmm->pi[i] = individual[num].number[phmm->N*(phmm->N+phmm->M)+i-1]; 
	}


	//////*****************************************************************
    double probf,probb;
	
	double **alpa,**beta;
	alpa = (double **) dmatrix(1, HmmMaxT, 1, phmm->N);
	beta = (double **) dmatrix(1, HmmMaxT, 1, phmm->N);
	
	double temp=0;

	for(k=0;k<LearnN;k++)
	{
		Forward(phmm,T[k],O[k],alpa,&probf);
		Backward(phmm,T[k],O[k],beta,&probb);
		
		for(t=1;t<T[k]-1;t++)
			for(i=1;i<=phmm->N;i++)
				for(j=1;j<=phmm->M;j++)
				{
					temp+=alpa[T[k]][i]*phmm->A[i][j]*phmm->B[j][O[k][T[k]+1]]*beta[T[k]+1][j];
				}
	}


	individual[num].fitness=temp;

	free_dmatrix(alpa,1,HmmMaxT,1,phmm->N);
	free_dmatrix(beta,1,HmmMaxT,1,phmm->N);

}


void FV(int i)
{
	int j;
	double sum1=0,sum2=0,sum3=0;
	for(j=0;j<DIMENSION;j++)
	{
		if(j<pow(HmmN,2))
			sum1=sum1+individual[i].number[j];
		else if(j>(pow(HmmN,2)-1)&&j<HmmN*(HmmN+HmmM))
			sum2=sum2+individual[i].number[j];
		else
			sum3=sum3+individual[i].number[j];
	}
	for(j=0;j<DIMENSION;j++)
	{
		if(j<pow(HmmN,2))
			individual[i].number[j]=individual[i].number[j]/sum1;
		else if(j>(pow(HmmN,2)-1)&&j<HmmN*(HmmN+HmmM))
			individual[i].number[j]=individual[i].number[j]/sum2;
		else
			individual[i].number[j]=individual[i].number[j]/sum3;
	}
	
	
}