Commit c5128c13 authored by Markus Mirz's avatar Markus Mirz
Browse files

rm old syn gen models in Archive folder

parent 3e79319f
#include "SynchronGenerator.h"
using namespace DPsim;
SynchronGenerator::SynchronGenerator(std::string name, int node1, int node2, int node3,
SynchGenStateType stateType, Real nomPower, Real nomVolt, Real nomFreq, int poleNumber, Real nomFieldCur,
SynchGenParamType paramType, Real Rs, Real Ll, Real Lmd, Real Lmd0, Real Lmq, Real Lmq0,
Real Rfd, Real Llfd, Real Rkd, Real Llkd,
Real Rkq1, Real Llkq1, Real Rkq2, Real Llkq2,
Real inertia) {
this->mNode1 = node1 - 1;
this->mNode2 = node2 - 1;
this->mNode3 = node3 - 1;
mStateType = stateType;
mNomPower = nomPower;
mNomVolt = nomVolt;
mNomFreq = nomFreq;
mPoleNumber = poleNumber;
mNomFieldCur = nomFieldCur;
// base stator values
mBase_V_RMS = mNomVolt / sqrt(3);
mBase_v = mBase_V_RMS * sqrt(2);
mBase_I_RMS = mNomPower / (3 * mBase_V_RMS);
mBase_i = mBase_I_RMS * sqrt(2);
mBase_Z = mBase_v / mBase_i;
mBase_OmElec = 2 * DPS_PI * mNomFreq;
mBase_OmMech = mBase_OmElec / (mPoleNumber / 2);
mBase_L = mBase_Z / mBase_OmElec;
mBase_Psi = mBase_L * mBase_i;
mBase_T = mNomPower / mBase_OmMech;
if (paramType == SynchGenParamType::perUnit) {
// steady state per unit initial value
initWithPerUnitParam(Rs, Ll, Lmd, Lmd0, Lmq, Lmq0, Rfd, Llfd, Rkd, Llkd, Rkq1, Llkq1, Rkq2, Llkq2, inertia);
}
}
void SynchronGenerator::initWithPerUnitParam(
Real Rs, Real Ll, Real Lmd, Real Lmd0, Real Lmq, Real Lmq0,
Real Rfd, Real Llfd, Real Rkd, Real Llkd,
Real Rkq1, Real Llkq1, Real Rkq2, Real Llkq2,
Real H) {
// base rotor values
mBase_ifd = Lmd * mNomFieldCur;
mBase_vfd = mNomPower / mBase_ifd;
mBase_Zfd = mBase_vfd / mBase_ifd;
mBase_Lfd = mBase_Zfd / mBase_OmElec;
if (mStateType == SynchGenStateType::perUnit) {
mRs = Rs;
mLl = Ll;
mLmd = Lmd;
mLmd0 = Lmd0;
mLmq = Lmq;
mLmq0 = Lmq0;
mRfd = Rfd;
mLlfd = Llfd;
mRkd = Rkd;
mLlkd = Llkd;
mRkq1 = Rkq1;
mLlkq1 = Llkq1;
mRkq2 = Rkq2;
mLlkq2 = Llkq2;
mH = H;
// Additional inductances according to Krause
mLaq = 1 / (1 / mLmq + 1 / mLl + 1 / mLlkq1 + 1 / mLlkq2);
mLad = 1 / (1 / mLmd + 1 / mLl + 1 / mLlkd + 1 / mLlfd);
}
else if (mStateType == SynchGenStateType::statorReferred) {
mRs = Rs * mBase_Z;
mLl = Ll * mBase_L;
mLmd = Lmd * mBase_L;
mLmd0 = Lmd0 * mBase_L;
mLmq = Lmq * mBase_L;
mLmq0 = Lmq0 * mBase_L;
mRfd = Rfd * mBase_Z;
mLlfd = Llfd * mBase_L;
mRkd = Rkd * mBase_Z;
mLlkd = Llkd * mBase_L;
mRkq1 = Rkq1 * mBase_Z;
mLlkq1 = Llkq1 * mBase_L;
mRkq2 = Rkq2 * mBase_Z;
mLlkq2 = Llkq2 * mBase_L;
// Additional inductances according to Krause
mLaq = 1 / (1 / mLmq + 1 / mLl + 1 / mLlkq1 + 1 / mLlkq2) * mBase_L;
mLad = 1 / (1 / mLmd + 1 / mLl + 1 / mLlkd + 1 / mLlfd) * mBase_L;
}
}
void SynchronGenerator::init(Real om, Real dt,
Real initActivePower, Real initReactivePower, Real initTerminalVolt, Real initVoltAngle) {
// Create matrices for state space representation
mInductanceMat <<
mLl + mLmq, 0, 0, mLmq, mLmq, 0, 0,
0, mLl + mLmd, 0, 0, 0, mLmd, mLmd,
0, 0, mLl, 0, 0, 0, 0,
mLmq, 0, 0, mLlkq1 + mLmq, mLmq, 0, 0,
mLmq, 0, 0, mLmq, mLlkq2 + mLmq, 0, 0,
0, mLmd, 0, 0, 0, mLlfd + mLmd, mLmd,
0, mLmd, 0, 0, 0, mLmd, mLlkd + mLmd;
mResistanceMat <<
mRs, 0, 0, 0, 0, 0, 0,
0, mRs, 0, 0, 0, 0, 0,
0, 0, mRs, 0, 0, 0, 0,
0, 0, 0, mRkq1, 0, 0, 0,
0, 0, 0, 0, mRkq2, 0, 0,
0, 0, 0, 0, 0, mRfd, 0,
0, 0, 0, 0, 0, 0, mRkd;
mOmegaFluxMat <<
0, 1, 0, 0, 0, 0, 0,
-1, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0;
mReverseCurrents <<
-1, 0, 0, 0, 0, 0, 0,
0, -1, 0, 0, 0, 0, 0,
0, 0, -1, 0, 0, 0, 0,
0, 0, 0, 1, 0, 0, 0,
0, 0, 0, 0, 1, 0, 0,
0, 0, 0, 0, 0, 1, 0,
0, 0, 0, 0, 0, 0, 1;
mReactanceMat = mInductanceMat.inverse();
if (mStateType == SynchGenStateType::perUnit) {
// steady state per unit initial value
initStatesInPerUnit(initActivePower, initReactivePower, initTerminalVolt, initVoltAngle);
}
else if (mStateType == SynchGenStateType::statorReferred) {
// steady state stator referred initial value
//InitStatesInStatorRefFrame(initActivePower, initReactivePower, initTerminalVolt, initVoltAngle);
}
mDq0Voltages(0, 0) = mVoltages(0, 0);
mDq0Voltages(1, 0) = mVoltages(1, 0);
mDq0Voltages(2, 0) = mVoltages(2, 0);
mDq0Voltages = mDq0Voltages * mBase_v;
mAbcsVoltages = dq0ToAbcTransform(mThetaMech, mDq0Voltages);
mDq0Currents(0, 0) = mCurrents(0, 0);
mDq0Currents(1, 0) = mCurrents(1, 0);
mDq0Currents(2, 0) = mCurrents(2, 0);
mDq0Currents = mDq0Currents * mBase_i;
mAbcsCurrents = dq0ToAbcTransform(mThetaMech, mDq0Currents);
}
void SynchronGenerator::initStatesInPerUnit(Real initActivePower, Real initReactivePower,
Real initTerminalVolt, Real initVoltAngle) {
double init_P = initActivePower / mNomPower;
double init_Q = initReactivePower / mNomPower;
double init_S = sqrt(pow(init_P, 2.) + pow(init_Q, 2.));
double init_vt = initTerminalVolt / mBase_v;
double init_it = init_S / init_vt;
// power factor
double init_pf = acos(init_P / init_S);
// load angle
double init_delta = atan(((mLmq + mLl) * init_it * cos(init_pf) - mRs * init_it * sin(init_pf)) /
(init_vt + mRs * init_it * cos(init_pf) + (mLmq + mLl) * init_it * sin(init_pf)));
double init_delta_deg = init_delta / DPS_PI * 180;
// dq stator voltages and currents
double init_vd = init_vt * sin(init_delta);
double init_vq = init_vt * cos(init_delta);
double init_id = init_it * sin(init_delta + init_pf);
double init_iq = init_it * cos(init_delta + init_pf);
// rotor voltage and current
double init_ifd = (init_vq + mRs * init_iq + (mLmd + mLl) * init_id) / mLmd;
double init_vfd = mRfd * init_ifd;
// flux linkages
double init_psid = init_vq + mRs * init_iq;
double init_psiq = -init_vd - mRs * init_id;
double init_psifd = (mLmd + mLlfd) * init_ifd - mLmd * init_id;
double init_psid1 = mLmd * (init_ifd - init_id);
double init_psiq1 = -mLmq * init_iq;
double init_psiq2 = -mLmq * init_iq;
// rotor mechanical variables
double init_Te = init_P + mRs * pow(init_it, 2.);
mOmMech = 1;
mVoltages(0, 0) = init_vq;
mVoltages(1, 0) = init_vd;
mVoltages(2, 0) = 0;
mVoltages(3, 0) = 0;
mVoltages(4, 0) = 0;
mVoltages(5, 0) = init_vfd;
mVoltages(6, 0) = 0;
mCurrents(0, 0) = init_iq;
mCurrents(1, 0) = init_id;
mCurrents(2, 0) = 0;
mCurrents(3, 0) = 0;
mCurrents(4, 0) = 0;
mCurrents(5, 0) = init_ifd;
mCurrents(6, 0) = 0;
mFluxes(0, 0) = init_psiq;
mFluxes(1, 0) = init_psid;
mFluxes(2, 0) = 0;
mFluxes(3, 0) = init_psiq1;
mFluxes(4, 0) = init_psiq2;
mFluxes(5, 0) = init_psifd;
mFluxes(6, 0) = init_psid1;
// Initialize mechanical angle
mThetaMech = initVoltAngle + init_delta - PI / 2.;
}
void SynchronGenerator::step(SystemModel& system, Real fieldVoltage, Real mechPower) {
if (mStateType == SynchGenStateType::perUnit) {
stepInPerUnit(system.getOmega(), system.getTimeStep(), fieldVoltage, mechPower);
}
else if (mStateType == SynchGenStateType::statorReferred) {
//StepInStatorRefFrame(om, dt, t, fieldVoltage, mechPower);
}
// Update current source accordingly
if (mNode1 >= 0) {
system.addCompToRightSideVector(mNode1, mAbcsCurrents(0, 0), mAbcsCurrents(3, 0));
}
if (mNode2 >= 0) {
system.addCompToRightSideVector(mNode2, mAbcsCurrents(1, 0), mAbcsCurrents(4, 0));
}
if (mNode3 >= 0) {
system.addCompToRightSideVector(mNode3, mAbcsCurrents(2, 0), mAbcsCurrents(5, 0));
}
}
void SynchronGenerator::stepInPerUnit(Real om, Real dt, Real fieldVoltage, Real mechPower) {
// retrieve voltages
mAbcsVoltages = (1 / mBase_v) * mAbcsVoltages;
mAbcsCurrents = (1 / mBase_i) * mAbcsCurrents;
// mVoltages(5, 0) = fieldVoltage / mBase_v;
// TODO calculate effect of changed field voltage
// dq-transform of interface voltage
mDq0Voltages = abcToDq0Transform(mThetaMech, mAbcsVoltages);
mVoltages(0, 0) = mDq0Voltages(0, 0);
mVoltages(1, 0) = mDq0Voltages(1, 0);
mVoltages(2, 0) = mDq0Voltages(2, 0);
// calculate mechanical states
mMechPower = mechPower / mNomPower;
mMechTorque = mMechPower / mOmMech;
mElecTorque = (mFluxes(1, 0)*mCurrents(0, 0) - mFluxes(0, 0)*mCurrents(1, 0));
// Euler step forward
mOmMech = mOmMech + dt * (1 / (2 * mH) * (mMechTorque - mElecTorque));
Matrix currents = mReverseCurrents * mReactanceMat * mFluxes;
DPSMatrix dtFluxes = mVoltages - mResistanceMat * currents - mOmMech * mOmegaFluxMat * mFluxes;
mFluxes = mFluxes + dt * mBase_OmElec * dtFluxes;
mCurrents = mReverseCurrents * mReactanceMat * mFluxes;
// inverse dq-transform
mDq0Currents(0, 0) = mCurrents(0, 0);
mDq0Currents(1, 0) = mCurrents(1, 0);
mDq0Currents(2, 0) = mCurrents(2, 0);
mAbcsCurrents = dq0ToAbcTransform(mThetaMech, mDq0Currents);
mAbcsCurrents = mBase_i * mAbcsCurrents;
// Update mechanical rotor angle with respect to electrical angle
mThetaMech = mThetaMech + dt * ((mOmMech - 1) * mBase_OmMech);
}
void SynchronGenerator::postStep(SystemModel& system) {
if (mNode1 >= 0) {
mAbcsVoltages(0, 0) = system.getRealFromLeftSideVector(mNode1);
mAbcsVoltages(3, 0) = system.getImagFromLeftSideVector(mNode1);
}
else {
mAbcsVoltages(0, 0) = 0;
mAbcsVoltages(3, 0) = 0;
}
if (mNode2 >= 0) {
mAbcsVoltages(1, 0) = system.getRealFromLeftSideVector(mNode2);
mAbcsVoltages(4, 0) = system.getImagFromLeftSideVector(mNode2);
}
else {
mAbcsVoltages(1, 0) = 0;
mAbcsVoltages(4, 0) = 0;
}
if (mNode3 >= 0) {
mAbcsVoltages(2, 0) = system.getRealFromLeftSideVector(mNode3);
mAbcsVoltages(5, 0) = system.getImagFromLeftSideVector(mNode3);
}
else {
mAbcsVoltages(2, 0) = 0;
mAbcsVoltages(5, 0) = 0;
}
}
DPSMatrix SynchronGenerator::abcToDq0Transform(Real theta, DPSMatrix& in) {
// Balanced case
Complex alpha(cos(2. / 3. * PI), sin(2. / 3. * PI));
Complex thetaCompInv(cos(-theta), sin(-theta));
MatrixComp AbcToPnz(3, 3);
AbcToPnz <<
1, 1, 1,
1, alpha, pow(alpha, 2),
1, pow(alpha, 2), alpha;
AbcToPnz = (1. / 3.) * AbcToPnz;
MatrixComp abcVector(3, 1);
abcVector <<
Complex(in(0, 0), in(3, 0)),
Complex(in(1, 0), in(4, 0)),
Complex(in(2, 0), in(5, 0));
MatrixComp pnzVector(3, 1);
pnzVector = AbcToPnz * abcVector * thetaCompInv;
DPSMatrix dq0Vector(3, 1);
dq0Vector <<
pnzVector(1, 0).imag(),
pnzVector(1, 0).real(),
0;
return dq0Vector;
}
DPSMatrix SynchronGenerator::dq0ToAbcTransform(Real theta, DPSMatrix& in) {
// Balanced case
Complex alpha(cos(2. / 3. * PI), sin(2. / 3. * PI));
Complex thetaComp(cos(theta), sin(theta));
MatrixComp PnzToAbc(3, 3);
PnzToAbc <<
1, 1, 1,
1, pow(alpha, 2), alpha,
1, alpha, pow(alpha, 2);
MatrixComp pnzVector(3, 1);
pnzVector <<
0,
Complex(in(1, 0), in(0, 0)),
Complex(0, 0);
MatrixComp abcCompVector(3, 1);
abcCompVector = PnzToAbc * pnzVector * thetaComp;
Matrix abcVector(6, 1);
abcVector <<
abcCompVector(0, 0).real(),
abcCompVector(1, 0).real(),
abcCompVector(2, 0).real(),
abcCompVector(0, 0).imag(),
abcCompVector(1, 0).imag(),
abcCompVector(2, 0).imag();
return abcVector;
}
#ifndef SYNCHRONGENERATOR_H
#define SYNCHRONGENERATOR_H
#include "BaseComponent.h"
#include "ComponentCommons.h"
namespace DPsim {
/// Synchronous generator model
/// If parInPerUnit is not set, the parameters have to be given with their respective stator or rotor
/// referred values. The calculation to per unit is performed in the initialization.
/// The case where parInPerUnit is not set will be implemented later.
/// parameter names include underscores and typical variables names found in literature instead of
/// descriptive names in order to shorten formulas and increase the readability
class SynchronGenerator : public BaseComponent {
protected:
// ### Machine parameters ###
/// nominal power Pn [VA]
Real mNomPower;
/// nominal voltage Vn [V] (RMS)
Real mNomVolt;
/// nominal frequency fn [Hz]
Real mNomFreq;
/// nominal field current Ifn [A]
Real mNomFieldCur;
/// stator resistance Rs[Ohm]
Real mRs;
/// leakage inductance Ll [H]
Real mLl;
/// d-axis mutual inductance Lmd [H]
Real mLmd;
/// unsaturated d-axis mutual inductance Lmd [H]
Real mLmd0;
/// q-axis mutual inductance Lmq [H]
Real mLmq;
/// unsaturated q-axis mutual inductance Lmq [H]
Real mLmq0;
/// field resistance Rfd [Ohm]
Real mRfd;
/// field leakage inductance Llfd [H]
Real mLlfd;
/// d-axis damper resistance Rkd [Ohm]
Real mRkd;
/// d-axis damper leakage inductance Llkd [H]
Real mLlkd;
/// q-axis damper resistance 1 Rkq1 [Ohm]
Real mRkq1;
/// q-axis damper leakage inductance 1 Llkq1 [H]
Real mLlkq1;
/// q-axis damper resistance 2 Rkq2 [Ohm]
Real mRkq2;
/// q-axis damper leakage inductance 2 Llkq2 [H]
Real mLlkq2;
/// q winding inductance
Real mLaq;
/// d winding inductance
Real mLad;
/// inertia J [kg*m^2]
Real mJ;
/// number of poles
int mPoleNumber;
/// inertia coefficient H
Real mH;
// ### Stator base values ###
/// specifies if the machine parameters are transformed to per unit
SynchGenStateType mStateType;
/// base stator voltage
Real mBase_v;
/// base stator voltage RMS
Real mBase_V_RMS;
/// base stator current
Real mBase_i;
/// base stator current RMS
Real mBase_I_RMS;
/// base stator impedance
Real mBase_Z;
/// base electrical angular frequency
Real mBase_OmElec;
/// base mechanical angular frequency
Real mBase_OmMech;
/// base stator inductance
Real mBase_L;
/// base torque
Real mBase_T;
/// base flux linkage
Real mBase_Psi;
/// base field current
Real mBase_ifd;
/// base field voltage
Real mBase_vfd;
/// base field impedance
Real mBase_Zfd;
/// base field inductance
Real mBase_Lfd;
// ### State variables ###
/// rotor speed omega_r
Real mOmMech;
/// theta
Real mThetaMech;
/// mechanical Power Pm [W]
Real mMechPower;
/// mechanical torque
Real mMechTorque;
/// electrical torque
Real mElecTorque;
/// voltage vector q d 0 kq1 kq2 df kd
DPSMatrix mVoltages = DPSMatrix::Zero(7, 1);
/// flux linkage vector
DPSMatrix mFluxes = DPSMatrix::Zero(7, 1);
/// current vector
DPSMatrix mCurrents = DPSMatrix::Zero(7, 1);
/// interface voltage vector abcs
DPSMatrix mAbcsVoltages = DPSMatrix::Zero(6, 1);
/// interface current vector abcs
DPSMatrix mAbcsCurrents = DPSMatrix::Zero(6, 1);
/// interface voltage vector dq0
DPSMatrix mDq0Voltages = DPSMatrix::Zero(3, 1);
/// interface current vector dq0
DPSMatrix mDq0Currents = DPSMatrix::Zero(3, 1);
// ### Useful Matrices ###
/// inductance matrix
DPSMatrix mInductanceMat = DPSMatrix::Zero(7, 7);
/// resistance matrix
DPSMatrix mResistanceMat = DPSMatrix::Zero(7, 7);
/// reactance matrix
DPSMatrix mReactanceMat = DPSMatrix::Zero(7, 7);
/// omega - flux matrix
DPSMatrix mOmegaFluxMat = DPSMatrix::Zero(7, 7);
/// matrix for reversing stator current directions in calculations with respect to other currents
DPSMatrix mReverseCurrents = DPSMatrix::Zero(7, 7);
public:
SynchronGenerator() { };
/// Initializes the per unit or stator referred machine parameters with the machine parameters given in per unit or
/// stator referred parameters depending on the setting of parameter type.
/// The initialization mode depends on the setting of state type.
SynchronGenerator(std::string name, int node1, int node2, int node3,
SynchGenStateType stateType, Real nomPower, Real nomVolt, Real nomFreq, int poleNumber, Real nomFieldCur,
SynchGenParamType paramType, Real Rs, Real Ll, Real Lmd, Real Lmd0, Real Lmq, Real Lmq0,
Real Rfd, Real Llfd, Real Rkd, Real Llkd,
Real Rkq1, Real Llkq1, Real Rkq2, Real Llkq2,
Real inertia);
/// Initializes the per unit or stator referred machine parameters with the machine parameters given in per unit.
/// The initialization mode depends on the setting of state type.
void initWithPerUnitParam(
Real Rs, Real Ll, Real Lmd, Real Lmd0, Real Lmq, Real Lmq0,
Real Rfd, Real Llfd, Real Rkd, Real Llkd, Real Rkq1, Real Llkq1,
Real Rkq2, Real Llkq2,
Real H);
/// Initializes states in per unit or stator referred variables depending on the setting of the state type.
/// Function parameters have to be given in Real units.
void init(Real om, Real dt,
Real initActivePower, Real initReactivePower, Real initTerminalVolt, Real initVoltAngle);
/// Initializes states in per unit. All machine parameters are assumed to be in per unit.
/// Function parameters have to be given in Real units.
void initStatesInPerUnit(Real initActivePower, Real initReactivePower,
Real initTerminalVolt, Real initVoltAngle);
/// Performs an Euler forward step with the state space model of a synchronous generator
/// to calculate the flux and current from the voltage vector.
void step(SystemModel& system, Real fieldVoltage, Real mechPower);
/// Performs an Euler forward step with the state space model of a synchronous generator
/// to calculate the flux and current from the voltage vector in per unit.
void stepInPerUnit(Real om, Real dt, Real fieldVoltage, Real mechPower);
/// Retrieves calculated voltage from simulation for next step
void postStep(SystemModel& system);
/// Park transform as described in Krause
DPSMatrix abcToDq0Transform(Real theta, DPSMatrix& in);
/// Inverse Park transform as described in Krause
DPSMatrix dq0ToAbcTransform(Real theta, DPSMatrix& in);
DPSMatrix getVoltages() { return mVoltages; }
DPSMatrix getCurrents() { return mCurrents; }
DPSMatrix getFluxes() { return mFluxes; }
void init(Real om, Real dt) { }
void applySystemMatrixStamp(SystemModel& system) { }
void applyRightSideVectorStamp(SystemModel& system) { }
};
}
#endif
#include "SynchronGeneratorEMT.h"
using namespace DPsim;
SynchronGeneratorEMT::SynchronGeneratorEMT(std::string name, int node1, int node2, int node3,
SynchGenStateType stateType, Real nomPower, Real nomVolt, Real nomFreq, int poleNumber, Real nomFieldCur,
SynchGenParamType paramType, Real Rs, Real Ll, Real Lmd, Real Lmd0, Real Lmq, Real Lmq0,
Real Rfd, Real Llfd, Real Rkd, Real Llkd,
Real Rkq1, Real Llkq1, Real Rkq2, Real Llkq2,
Real inertia) {
this->mNode1 = node1 - 1;
this->mNode2 = node2 - 1;
this->mNode3 = node3 - 1;
mStateType = stateType;
mNomPower = nomPower;