pso.h

求解梯级水库群优化调度

#include <stdlib.h>
#include <stdio.h>
#include <time.h>
#include <string>
#include <iostream>
#include <windows.h>
#include <vector>
#include <sstream>
#include "curve.h"
using namespace std;
#define PLANTNUM 3
#define SECTIONNUM 12
#define PSONUM 500
#define DIMNUM 36
#define A 500000000
#define B 6000000
#define C 800000
#define D 1000000
#define L 6000000000000
#define CONST_DECI 259.2
double gbest[DIMNUM];//群最好记录
double gbestdischarge[DIMNUM];//群最好记录对应弃水
double gbestobject;//群体最好解适应度函数值
double ibest[PSONUM][DIMNUM];//个体最好记录
double ibestdischarge[PSONUM][DIMNUM];//个体最好记录对应弃水
double ibestobject[PSONUM];//个体最好解适应度函数值
double maxp[DIMNUM]={1240,994,899,1240,994,899,1240,994,899,1240,994,899,1240,994,899,1240,994,899,1240,994,899,1240,994,899,1240,994,899,1240,994,899,1240,994,899,1236,994,899};//最大约束
double minp[DIMNUM]={1162,988,882,1162,988,882,1162,988,882,1162,988,882,1162,988,882,1162,988,882,1162,988,882,1162,988,882,1162,988,882,1162,988,882,1162,988,882,1162,988,882};//最小约束
//double interval[DIMNUM]={1550,63,207,2370,60,165,2640,57,577.5,2250,120,137.5,1590,38,164.8,874,25,156.4,553,18,82,420,11,58.5,380,10,33,418,10,2.85,597,8,15,908,7,35};
double interval[DIMNUM]={1580,20,70,2380,30,190,2720,30,300,2250,20,230,1560,20,190,889,10,131,559,7,74,419,4,55,375,4,38,401,5,26,585,6,16,847,10,25};
// double maxp[DIMNUM]={1240,994,899,1240,994,899,1240,994,899,1240,994,899};//最大约束
// double minp[DIMNUM]={1162,988,882,1162,988,882,1162,988,882,1162,988,882};//最小约束
// double interval[DIMNUM]={1550,63,207,2370,60,165,2640,57,577.5,2250,120,137.5};
//double slevel[PLANTNUM]={1190,990,894};//调度期初水位
double slevel[PLANTNUM]={1190,990,890};//调度期初水位
double qmax[PLANTNUM]={2200,2500,2142.6};//发电流量最大值
double qmin[PLANTNUM]={0,0,0};//发电流量最小值
double smax[PLANTNUM]={15,2,3};//水位最大速度
double maxpower[PLANTNUM]={4200000,1550000,1350000};
double minpower[PLANTNUM]={0,0,0};
double sprice[PLANTNUM];
double evagehead[PLANTNUM]={216,89,66};
double coeffiet[PLANTNUM]={8.44,8.6,8.3};
double evageflow[PLANTNUM]={1210,1230,1330};
double price[DIMNUM];
double speed[PSONUM][DIMNUM];//粒子飞行速度
double secpow[SECTIONNUM]={1000000,3000000,3000000,3500000,4000000,3300000,3200000,3000000,2900000,2600000,2200000,2000000};
double decimalrandom(double l,double h);//取l,h之间小数
int getrandom(int l,int h);//取l,h之间整数
void init(double **pso,double **discharge);//初始值
double mini(double a,double b);//两者之中取较小
double getpower(int t,int i,double *pso,double *discharge);//获取发电量
int fitness(double **pso,double **discharge,double *psovalue,string str);//适应度函数
double iterate(double **pso,double **discharge);//迭代
void update(double **pso,double **discharge,double *psovalue);//更新个体最佳，群体最佳
void upspeed(double w,double c1,double c2,double **pso);//更新速度
void GetDataOfString(vector<int > &vec,string str,char compart);
void GetDataOfString(vector<double > &vec,string str,char compart);
void stored_value();
void set_price(string str,int k);
string getstring (int n );
double getpowerkw(int t,int i,double *pso,double *discharge);
void stored_value()
{
	int i,j;
	double evage_price[PLANTNUM];
	for(i=0;i<PLANTNUM;i++)
		evage_price[i]=0;
	for(j=0;j<DIMNUM;j++)
		evage_price[j%PLANTNUM]+=price[j];
	for(j=0;j<PLANTNUM;j++)
		evage_price[j]/=SECTIONNUM;
    for(i=0;i<PLANTNUM;i++)
		sprice[i]=evage_price[i]*coeffiet[i]*evagehead[i]/3600;
	for(j=0;j<PLANTNUM;j++)
	{	
		for(i=j+1;i<PLANTNUM;i++)
			sprice[j]+=sprice[i];
		sprice[j]*=10000;
	}
}
void init(double **pso,double **discharge)
{
//初始化粒子及约束条件
	int i,j,m,t;
	double hmin,hmax;
	double v1,v2,down;
	for(i=0;i<PSONUM;i++)
	{
		for(j=0;j<DIMNUM;j++)
		{
			*(*(pso+i)+j)=0;
			*(*(discharge+i)+j)=0;
		}
	}
	for(i=0;i<PSONUM;i++)
	{
		down=0;
		for(j=0;j<DIMNUM;j++)
		{
			m=j%PLANTNUM;//电站
			t=j/PLANTNUM;//时段
// 			if(j>=DIMNUM-PLANTNUM)
// 			{
// 				*(*(pso+i)+j)=slevel[m];
// 				*(*(discharge+i)+j)=interval[j]+*(*(discharge+i)+j-1)-(getvolume(slevel[m],m)-getvolume(*(*(pso+i)+j-1),m))/CONST_DECI;
// 				*(*(ibestdischarge+i)+j)=*(*(discharge+i)+j);
// 			}
// 			else
// 			{
			hmin=*(minp+j);
			//下面初始解均可行
			v1=t>0 ? getvolume(*(*(pso+i)+j-PLANTNUM),m) : getvolume(*(slevel+m),m);
			v2=m>0 ? v1+(*(interval+j)+down)*CONST_DECI : v1+*(interval+j)*CONST_DECI;
			if(v2>getvolume(*(maxp+j),m))v2=getvolume(*(maxp+j),m);
			hmax=mini(*(maxp+j),getlevel(v2,m));
			*(*(pso+i)+j)=decimalrandom(hmin,hmax);
			*(*(ibest+i)+j)=*(*(pso+i)+j);
			v2=getvolume(*(*(pso+i)+j),m);
			down=m>0 ? *(interval+j)+(v1-v2)/CONST_DECI+down : *(interval+j)+(v1-v2)/CONST_DECI;
			*(*(discharge+i)+j)=down;
			*(*(ibestdischarge+i)+j)=*(*(discharge+i)+j);
// 			}
		}
	}
}
int fitness(double **pso,double **discharge,double *psovalue,string str)
{
//计算适应度函数值
	int i,t,j;
	double powerobject,storedobject;
	double *pp,*pd;
	double temp,tempv,tempp;
	double h,q,p;
    double weight_p=1,weight_s=1;
	double test,test_p;
	char strtemp[260];
	vector<double > d_temp;
	GetPrivateProfileString("Option","weight_p",NULL,strtemp,80,str.c_str());
	GetDataOfString(d_temp,strtemp,',');
	weight_p=d_temp.at(0);
	GetPrivateProfileString("Option","weight_s",NULL,strtemp,80,str.c_str());
	GetDataOfString(d_temp,strtemp,',');
	weight_s=d_temp.at(0);
	for(i=0;i<PSONUM;i++)
	{
		pp=*(pso+i);
	    pd=*(discharge+i);
		powerobject=0;
		storedobject=0;
		h=0;
		q=0;
		p=0;
	//	test=0;
	//	test_p=0;
		for(t=0;t<SECTIONNUM;t++)
		{
			for(j=0;j<PLANTNUM;j++)
			{
				temp=getpower(t,j,pp,pd);
				tempp=getpowerkw(t,j,pp,pd);
				powerobject+=temp*price[PLANTNUM*t+j];
				if(*(*(pso+i)+PLANTNUM*t+j)>*(maxp+PLANTNUM*t+j)+0.005)
					h+=*(*(pso+i)+PLANTNUM*t+j)-*(maxp+PLANTNUM*t+j)-0.005;
				else if(*(*(pso+i)+PLANTNUM*t+j)<*(minp+PLANTNUM*t+j)-0.005)
					h+=*(minp+PLANTNUM*t+j)-*(*(pso+i)+PLANTNUM*t+j)-0.005;
				if(*(*(discharge+i)+PLANTNUM*t+j)>*(qmax+j)+200)
					q+=*(*(discharge+i)+PLANTNUM*t+j)-*(qmax+j)-200;
				else if(*(*(discharge+i)+PLANTNUM*t+j)<*(qmin+j))
					q+=*(qmin+j)-*(*(discharge+i)+PLANTNUM*t+j);
				if(tempp>*(maxpower+j)+100)p+=tempp-*(maxpower+j)-100;
				else if(tempp<*(minpower+j))p+=*(minpower+j)-tempp;
			//	test+=tempp;
			}	
		//	test_p+=fabs(test-secpow[t]);			
		}
		for(j=0;j<PLANTNUM;j++)
		{
			tempv=getvolume(*(*(pso+i)+(SECTIONNUM-1)*(PLANTNUM)+j),j)-getvolume(*(slevel+j),j);
			storedobject+=tempv*sprice[j];
		}
	//	storedobject=0;
		*(psovalue+i)=L+weight_p*powerobject+weight_s*storedobject-A*h-B*q-C*p-D*test_p;
		*(psovalue+i)=L+weight_p*powerobject+weight_s*storedobject-A*h-B*q-C*p;
	}
// 	FILE *fp;
// 
// 	if((fp=fopen("E:\\epar\\PSO\\PSO\\test.ini","a"))==NULL)
// 	{
// 		printf("can not open the test.ini file");
// 		exit(0);
// 	}
// 	for(i=0;i<100;i++)
// 		fprintf(fp,"*");
// 	fprintf(fp,"\r\n");
// 	for(i=0;i<100;i++)
// 		fprintf(fp,"*");
// 	fprintf(fp,"\r\n");
// 	for(i=0;i<PSONUM;i++)
// 		fprintf(fp,"%f\r\n",*(psovalue+i));
	return 1;
}
void update(double **pso,double **discharge,double *psovalue)
{
//更新粒子最好解与群体最好解
	int i,j;
	for(i=0;i<PSONUM;i++)
	{
		if(*(psovalue+i)>gbestobject)
		{
			//更新群体最好记录
			gbestobject=*(psovalue+i);
			for(j=0;j<DIMNUM;j++)
			{
				*(gbest+j)=*(*(pso+i)+j);
				*(gbestdischarge+j)=*(*(discharge+i)+j);
			}
		}
		else if(*(psovalue+i)>*(ibestobject+i))
		{
			*(ibestobject+i)=*(psovalue+i);
			for(j=0;j<DIMNUM;j++)
			{
				*(*(ibest+i)+j)=*(*(pso+i)+j);
				*(*(ibestdischarge+i)+j)=*(*(discharge+i)+j);
			}
		}
	}
}
void upspeed(double w,double c1,double c2,double **pso)
{
	int i,j;
	double v1,v2,ib,gb;
	int m;
	for(i=0;i<PSONUM;i++)
	{
		for(j=0;j<DIMNUM;j++)
		{
// 			if(j>=DIMNUM-PLANTNUM)
// 				*(*(speed+i)+j)=0;
// 			else
// 			{
			v1=*(*(speed+i)+j);
			ib=*(*(ibest+i)+j);
			gb=*(gbest+j);
			v2=w*v1+c1*decimalrandom(0,1)*(ib-*(*(pso+i)+j))+c2*decimalrandom(0,1)*(gb-*(*(pso+i)+j));
			m=j%PLANTNUM;
			if(v2>*(smax+m))*(*(speed+i)+j)=*(smax+m);
			else if(v2<(-1)**(smax+m))*(*(speed+i)+j)=-1**(smax+m);
			else *(*(speed+i)+j)=v2;		
// 			}
		}
	}
}
double iterate(double **pso,double **discharge)
{
//迭代
	int i,j,t,m;
	double v1,v2;
	for(i=0;i<PSONUM;i++)
	{
		for(j=0;j<DIMNUM;j++)
		{
			t=j/PLANTNUM;
			m=j%PLANTNUM;
			*(*(pso+i)+j)+=*(*(speed+i)+j);
			v2=getvolume(*(*(pso+i)+j),m);
			if(t==0)
			{
				v1=getvolume(*(slevel+m),m);
				if(m==0)
					*(*(discharge+i)+j)=(v1-v2)/259.2+*(interval+j);
				else
					*(*(discharge+i)+j)=(v1-v2)/259.2+*(interval+j)+*(*(discharge+i)+j-1);
			}
			else
			{
				v1=getvolume(*(*(pso+i)+j-PLANTNUM),m);
				if(m==0)
					*(*(discharge+i)+j)=(v1-v2)/259.2+*(interval+j);
				else
					*(*(discharge+i)+j)=(v1-v2)/259.2+*(interval+j)+*(*(discharge+i)+j-1);
			}
		}
	}
// 	FILE *fp;
// 
// 	if((fp=fopen("E:\\epar\\PSO\\PSO\\level.ini","a"))==NULL)
// 	{
// 		printf("can not open the level.ini file");
// 		exit(0);
// 	}
// 	fprintf(fp,"\r\n");
// 	fprintf(fp,"\r\n");
// 	for(i=0;i<DIMNUM;i++)
// 		fprintf(fp,"%f\r\n",*(*(pso+0)+i));
 	return 0;
}
int getrandom(int l,int h)
{
//	调用函数之前请先调用一次srand((unsigned)time(NULL));
	return l+rand()%(h-l+1);
}
double decimalrandom(double l,double h)
{
// 	调用函数之前请先调用一次srand((unsigned)time(NULL));
	return (double)rand()/RAND_MAX*(h-l)+l;	
}
double mini(double a,double b)
{
	return a<b?a:b;
}
double getpowerkw(int t,int i,double *pso,double *discharge)
{
	double powerkw;
	double q,h,v,tail,down;
	if(t<1)
		v=(getvolume(*(slevel+i),i)+getvolume(*(pso+i),i))/2;
	else
		v=(getvolume(*(pso+PLANTNUM*(t-1)+i),i)+getvolume(*(pso+PLANTNUM*t+i),i))/2;
	down=*(discharge+PLANTNUM*t+i);
	tail=gettail(down,i);
	h=getlevel(v,i)-tail;
	q=down>qmax[i]?qmax[i]:down;	
    powerkw=coeffiet[i]*q*h;
	return powerkw;
}
double getpower(int t,int i,double *pso,double *discharge)//i是电站
{
	double power;
	double q,h,v,tail,down;
	if(t<1)
		v=(getvolume(*(slevel+i),i)+getvolume(*(pso+i),i))/2;
	else
		v=(getvolume(*(pso+PLANTNUM*(t-1)+i),i)+getvolume(*(pso+PLANTNUM*t+i),i))/2;
	down=*(discharge+PLANTNUM*t+i);
	tail=gettail(down,i);
	h=getlevel(v,i)-tail;
	q=down>qmax[i]?qmax[i]:down;	
    power=coeffiet[i]*q*h*720;
    return power;
}
void set_price(string str,int k)
{
	char strtemp[260];
	int j;
	vector<double > d_temp;
	double plant_p[PLANTNUM];
	string str_name="P";
	str_name+=getstring(k);
	GetPrivateProfileString("Price",str_name.c_str(),NULL,strtemp,80,str.c_str());
	GetDataOfString(d_temp,strtemp,',');
	for(j=0;j<PLANTNUM;j++)*(plant_p+j)=d_temp.at(j);
	for(j=0;j<DIMNUM;j++)price[j]=*(plant_p+j%PLANTNUM);
}
void GetDataOfString(vector<int > &vec,string str,char compart)
{
	int id;
	char temp[260];
    int start=0;
	int i,j,k;
	vec.clear();
	for(i=0;str[i]!='\0';i++)
	{
		if(str[i]==compart)
		{
			for(j=start,k=0;j<i;j++,k++)
				temp[k]=str[j];
			temp[k]='\0';
			id=atoi(temp);
			vec.push_back(id);
	    	start=i+1;
		}
	}
	for(j=start,k=0;j<i;j++,k++)
		temp[k]=str[j];
	temp[k]='\0';
	id=atoi(temp);
	vec.push_back(id);
}

void GetDataOfString(vector<double > &vec,string str,char compart)
{
	double d;
	char temp[260];
    int start=0;
	int i,j,k;
	vec.clear();
	for(i=0;str[i]!='\0';i++)
	{
		if(str[i]==compart)
		{
			for(j=start,k=0;j<i;j++,k++)
				temp[k]=str[j];
			temp[k]='\0';
			d=atof(temp);
			vec.push_back(d);
	    	start=i+1;
		}
	}
	for(j=start,k=0;j<i;j++,k++)
		temp[k]=str[j];
	temp[k]='\0';
	d=atof(temp);
	vec.push_back(d);
}
string getstring (int n )
{
	std::stringstream newstr;
    newstr<<n;
    return newstr.str();
}