exercise.cpp

这是一个龙格库塔算法计算的源代码

#include "stdafx.h"
#include "stdio.h"
#include "math.h"
#define PI 3.1415926

#if !defined(_CRunge_Kutta_H_)
#define _CRunge_Kutta_H_

class CRunge_Kutta
{
//attributions
private:
	double voltage;
	double Ids0;
	double Iqs0;
	double Idr0;
	double Iqr0;

	double Xss;
	double Xrr;
	double Xm;
	double R1;
	double R2;


	double Radian0;
	double wSpeed0;
	double wInertia;
	double wClose0;
	double Kd;
	double Time;

	double Tem;
	double Tl;

	double dT;
public:
//operation
private:
public:
	void SetParameter(double X1l,double X2l,double Xm0,
		double R1l,double R2l);
	void SetInitValue(double ids,double iqs,double idr,
		double iqr,double wSpeed,double Radian);//设置迭代初值 
	void SetJ(double value);//设置转动惯量
	void SetVoltage(double value);//设置电压
	void SetRadian(double value);//设置合闸角

	void SetTime(double value);//设置计算时间
	void SetdT(double value);//设置步长
	void SetKd(double value);//设置阻尼系数

	void Runge();//计算

	CRunge_Kutta();
	~CRunge_Kutta();
};
#endif

CRunge_Kutta::CRunge_Kutta()
{
	dT=0.001;
	Idr0=0;
	Iqr0=0;
	Ids0=0;
	Iqs0=0;
	Radian0=0;
	wClose0=0;
	Kd=0;
}
CRunge_Kutta::~CRunge_Kutta()
{
}

void CRunge_Kutta::SetdT(double value)
{
	dT=value;
}
void CRunge_Kutta::SetJ(double value)
{
	wInertia=4*value*2*PI*PI*PI*50*1500*1500/(60*60*2.2);
}
void CRunge_Kutta::SetKd(double value)
{
	Kd=value;
}
void CRunge_Kutta::SetTime(double value)
{
	Time=value;
}
void CRunge_Kutta::SetVoltage(double value)
{
	voltage=value;
}
void CRunge_Kutta::SetRadian(double value)
{
	wClose0=value;
}
void CRunge_Kutta::SetInitValue(double ids,double iqs,
								double idr,double iqr,
								double wSpeed,double Radian)
{
	Ids0=ids;
	Iqs0=iqs;
	Idr0=idr;
	Iqr0=iqr;
	wSpeed0=wSpeed;
	Radian0=Radian;
}

void CRunge_Kutta::SetParameter(double X1l,double X2l,
								double Xm0,double R1l,double R2l)
{
	Xm=Xm0;
	R1=R1l;
	R2=R2l;
	Xss=Xm0+X1l;
	Xrr=Xm0+X2l;
}

void CRunge_Kutta::Runge()
{
	int n=int(Time/dT);


	dT=100*PI*dT;//求时间标幺
	FILE* fp=fopen("Result.dat","w");
	for(int i=0;i<=n;++i)
	{
		double K1,K2,K3,K4;
		double Ids1,Iqs1,Idr1,Iqr1,wSpeed1,Radian1;
		//计算电流
		//Ids
		K1=(Xrr*(voltage*cos(i*dT+wClose0-Radian0)+wSpeed0*Xss*
			Iqs0+wSpeed0*Xm*Iqr0-R1*Ids0)+Xm*R2*Idr0)/(
			Xss*Xrr-Xm*Xm);
		K2=(Xrr*(voltage*cos(i*dT+dT/2+wClose0-Radian0)+wSpeed0*Xss*
			Iqs0+wSpeed0*Xm*Iqr0-R1*(Ids0+dT/2*K1))+
			Xm*R2*Idr0)/(Xss*Xrr-Xm*Xm);
		K3=(Xrr*(voltage*cos(i*dT+dT/2+wClose0-Radian0)+wSpeed0*Xss*
			Iqs0+wSpeed0*Xm*Iqr0-R1*(Ids0+dT/2*K2))+
			Xm*R2*Idr0)/(Xss*Xrr-Xm*Xm);
		K4=(Xrr*(voltage*cos(i*dT+dT+wClose0-Radian0)+wSpeed0*Xss*
			Iqs0+wSpeed0*Xm*Iqr0-R1*(Ids0+dT*K3))+
			Xm*R2*Idr0)/(Xss*Xrr-Xm*Xm);
		Ids1=dT*(K1+2*K2+2*K3+K4)/6.0;
		//Iqs
		K1=(Xrr*(voltage*sin(i*dT+wClose0-Radian0)-wSpeed0*Xss*
			Ids0-wSpeed0*Xm*Idr0-R1*Iqs0)+Xm*R2*Iqr0)/(
			Xss*Xrr-Xm*Xm);
		K2=(Xrr*(voltage*sin(i*dT+dT/2+wClose0-Radian0)-wSpeed0*Xss*
			Ids0-wSpeed0*Xm*Idr0-R1*(Iqs0+dT/2*K1))+
			Xm*R2*Iqr0)/(Xss*Xrr-Xm*Xm);
		K3=(Xrr*(voltage*sin(i*dT+dT/2+wClose0-Radian0)-wSpeed0*Xss*
			Ids0-wSpeed0*Xm*Idr0-R1*(Iqs0+dT/2*K2))+
			Xm*R2*Iqr0)/(Xss*Xrr-Xm*Xm);
		K4=(Xrr*(voltage*sin(i*dT+dT+wClose0-Radian0)-wSpeed0*Xss*
			Ids0-wSpeed0*Xm*Idr0-R1*(Iqs0+dT*K3))+
			Xm*R2*Iqr0)/(Xss*Xrr-Xm*Xm);
		Iqs1=dT*(K1+2*K2+2*K3+K4)/6.0;
		//Idr
		K1=-(Xm*(voltage*cos(i*dT+wClose0-Radian0)+wSpeed0*Xss*
			Iqs0+wSpeed0*Xm*Iqr0-R1*Ids0)+Xss*R2*Idr0)/(
			Xss*Xrr-Xm*Xm);
		K2=-(Xm*(voltage*cos(i*dT+dT/2+wClose0-Radian0)+wSpeed0*Xss*
			Iqs0+wSpeed0*Xm*Iqr0-R1*Ids0)+
			Xss*R2*(Idr0+dT/2*K1))/(Xss*Xrr-Xm*Xm);
		K3=-(Xm*(voltage*cos(i*dT+dT/2+wClose0-Radian0)+wSpeed0*Xss*
			Iqs0+wSpeed0*Xm*Iqr0-R1*Ids0)+
			Xss*R2*(Idr0+dT/2*K2))/(Xss*Xrr-Xm*Xm);
		K4=-(Xm*(voltage*cos(i*dT+dT+wClose0-Radian0)+wSpeed0*Xss*
			Iqs0+wSpeed0*Xm*Iqr0-R1*Ids0)+
			Xss*R2*(Idr0+dT*K3))/(Xss*Xrr-Xm*Xm);
		Idr1=dT*(K1+2*K2+2*K3+K4)/6.0;
		//Iqr
		K1=-(Xm*(voltage*sin(i*dT+wClose0-Radian0)-wSpeed0*Xss*
			Ids0-wSpeed0*Xm*Idr0-R1*Iqs0)+Xss*R2*Iqr0)/(
			Xss*Xrr-Xm*Xm);
		K2=-(Xm*(voltage*sin(i*dT+dT/2+wClose0-Radian0)-wSpeed0*Xss*
			Ids0-wSpeed0*Xm*Idr0-R1*Iqs0)+
			Xss*R2*(Iqr0+dT/2*K1))/(Xss*Xrr-Xm*Xm);
		K3=-(Xm*(voltage*sin(i*dT+dT/2+wClose0-Radian0)-wSpeed0*Xss*
			Ids0-wSpeed0*Xm*Idr0-R1*Iqs0)+
			Xss*R2*(Iqr0+dT/2*K2))/(Xss*Xrr-Xm*Xm);
		K4=-(Xm*(voltage*sin(i*dT+dT+wClose0-Radian0)-wSpeed0*Xss*
			Ids0-wSpeed0*Xm*Idr0-R1*Iqs0)+
			Xss*R2*(Iqr0+dT*K3))/(Xss*Xrr-Xm*Xm);
		Iqr1=dT*(K1+2*K2+2*K3+K4)/6.0;

		//计算电磁转矩
		Tem=Xm*(Idr0*Iqs0-Iqr0*Ids0);
		Tl=0.007012*wSpeed0*wSpeed0;
		K1=(Tem-Tl-Kd*wSpeed0)/wInertia;
		Tl=0.007012*(wSpeed0+dT*K1/2)*(wSpeed0+dT*K1/2);
		K2=(Tem-Tl-Kd*(wSpeed0+dT*K1/2))/wInertia;
		Tl=0.007012*(wSpeed0+dT*K2/2)*(wSpeed0+dT*K2/2);
		K3=(Tem-Tl-Kd*(wSpeed0+dT*K2/2))/wInertia;
		Tl=0.007012*(wSpeed0+dT*K3)*(wSpeed0+dT*K3);
		K4=(Tem-Tl-Kd*(wSpeed0+dT*K3))/wInertia;

		wSpeed1=dT*(K1+2*K2+2*K3+K4)/6.0;

		//计算相角
		Radian1=dT*wSpeed0;

		double Ia=Ids0*cos(Radian0)-Iqs0*sin(Radian0);
		Ids0+=Ids1;
		Iqs0+=Iqs1;
		Idr0+=Idr1;
		Iqr0+=Iqr1;
		wSpeed0+=wSpeed1;
		Radian0+=Radian1;
		double value0=1430.0/1500.0;
		fprintf(fp,"%8.6f\n",wSpeed0-1.0);
//		fprintf(fp,"%d\t%8.6f\t%8.6f\t%8.6f\t%8.6f\t%8.6f\t%8.6f\t%8.6f\t%8.6f\n",
//			i+1,Ids0,Iqs0,Idr0,Iqr0,wSpeed0,Radian0,Tem,Ia);
	}
	fclose(fp);
}

int main(int argc, char* argv[])
{
	CRunge_Kutta* pRun=new CRunge_Kutta();
	FILE *fp=fopen("Data.dat","r");

	float value1,value2,value3,
		value4,value5,value6;
	fscanf(fp,"%f%f%f%f%f",&value1,&value2,
		&value3,&value4,&value5);
	pRun->SetParameter(value1,value2,
		value3,value4,value5);
	fscanf(fp,"%f%f%f%f%f%f",&value1,&value2,
		&value3,&value4,&value5,&value6);
	pRun->SetVoltage(value1);
	pRun->SetJ(value2);
	pRun->SetRadian(value3);
	pRun->SetKd(value4);
	pRun->SetTime(value5);
	pRun->SetdT(value6);
	fscanf(fp,"%f%f%f%f%f%f",&value1,&value2,
		&value3,&value4,&value5,&value6);
	pRun->SetInitValue(value1,value2,value3,
		value4,value5,value6);
	fclose(fp);

	pRun->Runge();
	return 0;
}