tabulist.cpp

求解梯级水库群优化调度

#ifndef _TABULIST_H_
#include "TabuList.h"
#ifndef _PLANT_H_
#include "Plant.h"
#endif

CTabuList::CTabuList()
{
	Init();
}
CTabuList::CTabuList(int iteratornum,int tabulength,int neighbornum)
{
	Init();
	IteratorNum=iteratornum;
	TabuLength=tabulength;
	NeighborNum=neighbornum;
}
CTabuList::~CTabuList()
{
	
}
void CTabuList::Init()
{
	IteratorNum=100;//初始化迭代次数
	TabuLength=20;//初始化禁忌长度
	NeighborNum=10;//初始化邻域解个数
	TabuList.clear();//禁忌表清空
	TabuListAspiration.clear();//禁忌表对应期望值清空
 	InitPlant();//初始化电站
	InitStep();//初始化步长
 	int num_t,num_p;//存放时段数和电站数
	num_t=CPlant::SectionNum;
	num_p=CPlant::PlantNum;
	deque<double > temp;
	for(int i=0;i<num_p;i++)
		temp.push_back(0);
	for(int j=0;j<num_t;j++)
	{//当前解
		CurrentWaterLevel.push_back(temp);
		CurrentPowerFlow.push_back(temp);
		CurrentDischarge.push_back(temp);
	}
    InitSolution();
	//当前最好解=当前解（初始解）
	BestSoFarWaterLevel=CurrentWaterLevel;
	BestSoFarPowerFlow=CurrentPowerFlow;
	BestSoFarDischarge=CurrentDischarge;
	StartWaterLevel=CurrentWaterLevel;
	StartPowerFlow=CurrentPowerFlow;
	StartDischarge=CurrentDischarge;
	CurrentObjectValue=ObjectValue(CurrentWaterLevel,CurrentPowerFlow,CurrentDischarge);//得到目标函数值
	BestSoFarObjectValue=ObjectValue(BestSoFarWaterLevel,BestSoFarPowerFlow,BestSoFarDischarge);
	StartObject=CurrentObjectValue;
	PushToTabuList(CurrentWaterLevel);//初始可行解放进紧急表
    TabuListAspiration.push_back(CurrentObjectValue);//初始解对应目标函数放进禁忌表对应期望值
	NeighborWaterLevel.clear();//清空邻域解及邻域解对应之目标函数值
	NeighborPowerFlow.clear();
	NeighborDischarge.clear();
	NeighborObjectValue.clear();
	cout<<"初始解结果："<<CurrentObjectValue<<endl;
}
void CTabuList::InitPlant()//初始化电站信息
{
	CPlant cpt1("XW","e:\\epar\\Rsos8.29\\test.ini");
// 	CPlant cpt2("MW","e:\\epar\\Rsos8.29\\test.ini");
// 	CPlant cpt3("JH","e:\\epar\\Rsos8.29\\test.ini");
	plants[0]=cpt1;
// 	plants[1]=cpt2;
// 	plants[2]=cpt3;
}

void CTabuList::InitStep()//初始化步长矢量
{
	int i;
	for(i=0;i<CPlant::PlantNum;i++)
	{
		double temp=plants[i].GetPlantStepLength();
		Step.push_back(temp);//初始化步长向量
		StepAttept.push_back(0);
		StepAccept.push_back(0);		
	}

}
void CTabuList::InitSolution()//初始可行解
{
	int i,j,flag;
	flag=1;
	for(i=0;i<CPlant::SectionNum;i++)
	{
		flag=1;
		for(j=0;j<CPlant::PlantNum;j++)
		{
			double sl,down,lastdown,pftemp,htemp,hmax,hmin,interval,pfmin,pfmax,spi;
			sl=0,down=0,lastdown=0,pftemp=0,htemp=0,hmax=0,hmin=0,interval=0,pfmin=0,pfmax=0,spi=0;
			sl=i>0?CurrentWaterLevel[i-1][j]:plants[j].GetStartLevel();
			lastdown=j>0?(CurrentPowerFlow[i][j-1]+CurrentDischarge[i][j-1]):0;
			interval=plants[j].GetSectionIntervalFlow(i);
			plants[j].GetMaxMinLevel(hmax,hmin,sl,lastdown,i);
// 			hmax=plants[j].GetMaxWaterLevel(i);
// 			hmin=plants[j].GetMinWaterLevel(i);
			pfmin=plants[j].GetMinPowerFlow(i);
			pfmax=plants[j].GetMaxPowerFlow(i);
			do 
			{
				htemp=getrandom(hmin,hmax);
			    down=plants[j].GetPlantDownFlow(sl,htemp,i,lastdown);
			    Sleep(10);
				hmax--;
				if(hmax<=hmin)
				{
					flag=0;
					break;
				}
			} while(down<pfmin);
			if(flag==0)
			{
				j--;
				i--;
				break;
			}
			pftemp=down>pfmax?pfmax:down;
			spi=down>pfmax?down-pfmax:0;
			CurrentWaterLevel[i][j]=htemp;
			CurrentPowerFlow[i][j]=down>pfmax?pfmax:down;
			CurrentDischarge[i][j]=down>pfmax?down-pfmax:0;
			cout.flags(ios::right);
			cout<<"期末水位："<<setw(8)<<htemp<<" 下泄流量："<<setw(8)<<down<<" 发电流量："<<setw(8)<<pftemp<<" 弃水流量："<<setw(8)<<spi<<endl;
		}
	}
	cout<<"the end "<<endl;
	for(i=0;i<CPlant::SectionNum;i++)
		for(j=0;j<CPlant::PlantNum;j++)
		{
			cout<<"期末水位："<<setw(8)<<CurrentWaterLevel[i][j]<<"发电流量："<<setw(8)<<CurrentPowerFlow[i][j];
		    cout<<"弃水流量："<<setw(8)<<CurrentDischarge[i][j]<<endl;
		}
}

void CTabuList::GetNeighborSolutions()
{
	//调整步长应该防御此
	NeighborWaterLevel.clear();
	NeighborPowerFlow.clear();
	NeighborDischarge.clear();
	NeighborObjectValue.clear();
	deque<deque<double > > cw,cp,cd;
	int n;
	int K,T;
	double htemp,hmax,hmin,dmax,dmin,pfmax,pfmin,pftemp,down,lastdown,interval,sl,t;
	for(n=0;n<NeighborNum;n++)
	{
	//	cout<<"the "<<n<<"-th neighbor"<<endl;	
		int flag=1;
		do 
		{
	        cw=CurrentWaterLevel;
	        cp=CurrentPowerFlow;
	        cd=CurrentDischarge;
			flag=1;
		//	cout<<"in the do while for the "<<n<<"-th neighbor"<<endl;	
			double r=0;
			K=getrandom(0,CPlant::PlantNum-1);
		//	Sleep(100);
			T=getrandom(0,CPlant::SectionNum-1);
	//		cout<<"时段是"<<T<<"电站是"<<K<<endl;////////////////////////////////////////////////
			for(int i=T;i<CPlant::SectionNum;i++)
			{
				for(int j=(i==T)?K:0;j<CPlant::PlantNum;j++)
				{		
					htemp=0,hmax=0,hmin=0,dmax=0,dmin=0,pfmax=0,pfmin=0,pftemp=0,down=0,lastdown=0,interval=0,sl=0;
					sl=i>0?cw[i-1][j]:plants[j].GetStartLevel();
					lastdown=j>0?cp[i][j-1]+cd[i][j-1]:0;
					plants[j].GetMaxMinLevel(hmax,hmin,sl,lastdown,i);
// 					hmax=plants[j].GetMaxWaterLevel(i);
// 					hmin=plants[j].GetMinWaterLevel(i);
					interval=plants[j].GetSectionIntervalFlow(i);
					dmin=plants[j].GetMinDownFlow(i);
					dmax=plants[j].GetMaxDownFLow(i);
					pfmax=plants[j].GetMaxPowerFlow(i);
					pfmin=plants[j].GetMinPowerFlow(i);
					t=plants[j].GetSectionTime();
					if(j==K && i==T)//对于随机产生的电站时段
					{
						r=getrandom(0.0,1.0)*pow(-1,getrandom(0,1));//调整系数（-1，1）
	//					cout<<"r 是："<<r<<endl;///////////////////////////////////////
						htemp=cw[i][j]+Step[j]*r;
						if(between(hmin,hmax,htemp))//若产生水位满足最大最小约束，则求下泄流量
						{
							double v1,v2;
							v1=plants[j].GetPlantVolume(sl);
							v2=plants[j].GetPlantVolume(htemp);			
							down=(v1-v2)*10000/t+interval+lastdown;
							if(!between(pfmin,pfmax,down))
								flag=flag&&0;//若下泄流量不满足发电最大最小约束，跳出循环
							else//若下泄流量满足最大最小约束，更新T时段K电站的水位，发电流量，泄水流量
							{
								cw[i][j]=htemp;
								cp[i][j]=down;
								cd[i][j]=0;
							}
						}
						else 
							flag=flag && 0;//若产生水位不满足最大最小约束，则跳出循环						
					}
					else//对与非随机产生的时段电站
					{
						double v1,v2;
						htemp=cw[i][j];
						v1=plants[j].GetPlantVolume(sl);
						v2=plants[j].GetPlantVolume(htemp);
						down=(v1-v2)*10000/t+interval+lastdown;
						if(down<=pfmax && down>=pfmin)
						{
							cw[i][j]=htemp;
							cp[i][j]=down;
							cd[i][j]=0;
						}
						else if(down>pfmax)
						{
						 	htemp=plants[j].GetPlantLevel(plants[j].GetPlantVolume(cw[i][j])+(down-pfmax)*t/10000);
							cp[i][j]=pfmax;
							if(htemp<=hmax)
							{
								cw[i][j]=htemp;
								cd[i][j]=0;
							}
							else
							{
								cw[i][j]=hmax;
								cd[i][j]+=(plants[j].GetPlantVolume(htemp)-plants[j].GetPlantVolume(hmax))*10000/t;
							//	cout<<"the down is "<<down<<endl;
							}
						}
						else if(down<pfmin)
						{
						//	cout<<"the down < pfmin down is :"<<down<<endl;//////////////
							htemp=plants[j].GetPlantLevel(plants[j].GetPlantVolume(cw[i][j])+(down-pfmin)*t/10000);
							cp[i][j]=pfmin;
							if(htemp>=hmin)
							{
								cw[i][j]=htemp;
								cd[i][j]=0;
							}
							else
							{
								cw[i][j]=hmin;
								cd[i][j]+=(plants[j].GetPlantVolume(htemp)-plants[j].GetPlantVolume(hmin))*10000/t;
							}
						}
					}
					if(!flag)
						break;
				}
				if(!flag)
					break;
			}
		//	Sleep(1000);
			StepAttept[K]++;
			for(i=0;i<CPlant::SectionNum;i++)
			{
				for(int j=0;j<CPlant::PlantNum;j++)
				{
					double v1,v2;
					t=plants[j].GetSectionTime();
					sl=i>0?cw[i-1][j]:plants[j].GetStartLevel();
					lastdown=j>0?cp[i][j-1]+cd[i][j-1]:0;
					interval=plants[j].GetSectionIntervalFlow(i);
					down=cp[i][j]+cd[i][j];
					v2=plants[j].GetPlantVolume(cw[i][j])*10000;
					v1=plants[j].GetPlantVolume(sl)*10000;
					if(fabs(v2-v1-(interval+lastdown-down)*t)>10000)
					{
					//	cout<<"dayu 10000"<<"the flag is :"<<flag<<endl;
						flag=flag && 0;
						break;
					}						
				}
				if(!flag)break;
			}
			//检测所有时段所有电站是否满足水量平衡,此步至关重要，若全部满足，则为可行解，否则为非可行解
			//若为可行解则accept加1
			////尝试次数加1
		//	if(flag)cout<<"tong guo"<<endl;
		} while(!flag);
//		cout<<"the "<<n<<"-th neighbor solution"<<endl;
		NeighborWaterLevel.push_back(cw);
		NeighborPowerFlow.push_back(cp);
		NeighborDischarge.push_back(cd);
//		cout<<"the dian zhan shi "<<K<<"the shi duan is "<<T<<endl;
		StepAccept[K]++;
// 		for(int s=0;s<CPlant::SectionNum;s++)
// 			for(int t=0;t<CPlant::PlantNum;t++)
// 				cout<<"水位 "<<cw[s][t]<<"发电 "<<cp[s][t]<<"弃水 "<<cd[s][t]<<endl;
	}
	
	for(int i=0;i<CPlant::PlantNum;i++)
		cout<<"尝试次数是 :"<<StepAttept[i]<<"接受次数是："<<StepAccept[i]<<endl;// 
	StepAdjust();

}

void CTabuList::PushToTabuList(deque<deque<double > > &tabuitem)
{
	if(TabuList.size()<TabuLength)
		TabuList.push_back(tabuitem);
	else
	{
		TabuList.pop_front();
		TabuList.push_back(tabuitem);
	}
}
void CTabuList::UpdateAspiration(double l)
{
	if(TabuListAspiration.size()<TabuLength)
		TabuListAspiration.push_back(l);
	else
	{
		TabuListAspiration.pop_front();
		TabuListAspiration.push_back(l);
	}
}
bool CTabuList::IsTabu(deque<deque<double > > &neighbor)
{
	bool b=false;
	for(int i=0;i<TabuList.size();i++)
	{
		if(neighbor==TabuList.at(i))
		{
			b=true;
			break;
		}
	}
	return b;
}
bool CTabuList::AlTest(int n)
{
	if(NeighborObjectValue.at(n)>BestSoFarObjectValue)
		return true;
	else
		return false;
}
void CTabuList::StepAdjust()
{
	for(int i=0;i<CPlant::PlantNum;i++)
	{
		if(StepAttept[i]!=0 && (double)StepAccept[i]/StepAttept[i]>P_MAX)
		{
			Step[i]=Step[i]*STEP_MAX;//是否有问题，有待考虑
		//	cout<<"the stepattempt i is :"<<StepAttept[i]<<endl;
		}
		else if(StepAttept[i]!=0 && (double)StepAccept[i]/StepAttept[i]<P_MIN)
		{
			Step[i]=Step[i]*STEP_MIN;
		//	cout<<"the stepattempt i is :"<<StepAttept[i]<<endl;
		}
		StepAttept[i]=0;
		StepAccept[i]=0;
	}
}

double CTabuList::ObjectValue(deque<deque<double > > &w,deque<deque<double > > &p,deque<deque<double > > &d)
{
	double powerprofit=0,storedvalue=0;
	double k,h,hf,hl,powerprice,t,v1,v2,sp,hdead;
	int i=0,j=0;
	for(i=0;i<CPlant::SectionNum;i++)
	{
		for(j=0;j<CPlant::PlantNum;j++)
		{
			k=plants[j].GetEvageCoeff();
			hf=plants[j].GetPlantTailLevel(p[i][j]+d[i][j]);
//			hl=plants[j].GetPlantLostLevel(p[i][j]);//目前水头损失函数未加，故采用平均水头损失
			hl=plants[j].GetEvageLevelLost();
            h=i>0?(w[i-1][j]+w[i][j])/2-hf-hl:(plants[j].GetStartLevel()+w[i][j])/2-hf-hl;
			powerprice=plants[j].GetSectionPowerPrice(i);
			t=plants[j].GetSectionTime();
			powerprofit+=powerprice*k*p[i][j]*h*t/3600;			
		}
	}
// 	cout<<"the power profit is : "<<powerprofit<<endl;
	for(i=0;i<CPlant::PlantNum;i++)
	{
		hdead=plants[i].GetDeadLevel();
		v1=plants[i].GetPlantVolume(hdead);
		v2=plants[i].GetPlantVolume(w[CPlant::SectionNum-1][i]);
		sp=plants[i].GetStoredPrice();
// 		cout<<v2-v1<<endl;
		storedvalue+=sp*(v2-v1)*10000;//存水应该仅仅到死水位
	}
// 	cout<<"the stored value is : "<<storedvalue<<endl;
	return powerprofit+storedvalue;
}


void CTabuList::OutputResult()
{
	int plant_num,section_num;
	string plant_name;
	plant_num=CPlant::PlantNum;
	section_num=CPlant::SectionNum;
	cout<<endl;
	cout<<endl;
	cout<<"本次仿真优化调度结果如下："<<endl;
	cout<<endl;
	cout<<"初始解目标是："<<StartObject<<endl;
	cout<<"目标函数值是："<<BestSoFarObjectValue<<endl;
	cout<<"初始及期间各电站期末水位、发电流量、弃水流量如下："<<endl;
    for(int i=0;i<plant_num;i++)
	{
		plant_name=plants[i].GetPlantName();
		cout<<endl;
		cout<<plant_name<<":"<<endl;
		for(int j=0;j<section_num;j++)
			cout<<BestSoFarWaterLevel[j][i]<<'\t'<<'\t'<<BestSoFarPowerFlow[j][i]<<'\t'<<'\t'<<BestSoFarDischarge[j][i]<<endl;
		cout<<endl;
	}
	cout<<"禁忌表长:"<<TabuList.size()<<endl;
}

void CTabuList::RunTabuSearch(int n)
{

}
void CTabuList::RunTabuSearch()
{ 
// 	srand((int)time(NULL));
// 	GetNeighborSolutions();
// 	GetNeighborSolutions();

	int i,j;
	//GetNeighborSolutions();
	int CurrentStepNum=0;
	double MaxObject,temp;
	int MaxNum;
	srand((int)time(NULL));
	for(i=0;i<IteratorNum;i++,CurrentStepNum++)
	{
		clock_t start,finish;
	
		start=clock();
		GetNeighborSolutions();
		finish=clock();
		cout<<"getneighbor time is :"<<(finish-start)/CLOCKS_PER_SEC<<endl;
		MaxObject=ObjectValue(NeighborWaterLevel.at(0),NeighborPowerFlow.at(0),NeighborDischarge.at(0));
		MaxNum=0;
		for(j=0;j<NeighborNum;j++)
		{
			if(IsTabu(NeighborWaterLevel.at(j)))
			{//是禁忌状态
				if(AlTest(j))
				{//满足特赦准则
					CurrentWaterLevel=NeighborWaterLevel.at(j);
				    CurrentPowerFlow=NeighborPowerFlow.at(j);
				    CurrentDischarge=NeighborDischarge.at(j);
				    BestSoFarWaterLevel=CurrentWaterLevel;
				    BestSoFarPowerFlow=CurrentPowerFlow;
				    BestSoFarDischarge=CurrentDischarge;
				    CurrentObjectValue=ObjectValue(CurrentWaterLevel,CurrentPowerFlow,CurrentDischarge);
				    BestSoFarObjectValue=CurrentObjectValue;
				}
			}
			else
			{//若不是禁忌状态
				temp=ObjectValue(NeighborWaterLevel.at(j),NeighborPowerFlow.at(j),NeighborDischarge.at(j));
				if(temp>CurrentObjectValue)
				{
					CurrentWaterLevel=NeighborWaterLevel.at(j);
					CurrentPowerFlow=NeighborPowerFlow.at(j);
					CurrentDischarge=NeighborDischarge.at(j);
					CurrentObjectValue=temp;
					if(CurrentObjectValue>BestSoFarObjectValue)
					{
						BestSoFarWaterLevel=CurrentWaterLevel;
						BestSoFarPowerFlow=CurrentPowerFlow;
						BestSoFarDischarge=CurrentDischarge;
						BestSoFarObjectValue=CurrentObjectValue;
						PushToTabuList(CurrentWaterLevel);
						UpdateAspiration(CurrentObjectValue);
					}
				}
			}
		}
		cout<<"第"<<i<<"次迭代目标函数值： "<<BestSoFarObjectValue<<endl;
	
	}
	OutputResult();
}

void CTabuList::SetTabuLength(int l)
{
	TabuLength=l;
}
void CTabuList::SetNeighborNum(int n)
{
	NeighborNum=n;
}

#endif