Commit 4088f495 authored by Viviane Sapucaia's avatar Viviane Sapucaia
Browse files

Fixed VBR Example

parent 08c76619
/** SynGenVBR Example
/** SynGenVBR Example
*
* @author Markus Mirz <mmirz@eonerc.rwth-aachen.de>
* @copyright 2017, Institute for Automation of Complex Power Systems, EONERC
......@@ -26,9 +26,6 @@ using namespace DPsim::Components::EMT;
int main(int argc, char* argv[])
{
// Define Object for saving data on a file
Logger vtLog("data_vt.csv"),
jLog("data_j.csv");
// Define machine parameters in per unit
Real nomPower = 555e6;
......@@ -51,35 +48,35 @@ int main(int argc, char* argv[])
Real Llkd = 0.1713;
Real Rkq1 = 0.0062;
Real Llkq1 = 0.7252;
//Real Rkq2 = 0.0237;
//Real Llkq2 = 0.125;
Real Rkq2 = 0;
Real Llkq2 = 0;
Real Rkq2 = 0.0237;
Real Llkq2 = 0.125;
//Real Rkq2 = 0;
//Real Llkq2 = 0;
// Declare circuit components
Components::Base::Ptr gen = SynchronGeneratorVBR::make("gen", 1, 2, 3,
Components::Base::Ptr gen = SynchronGeneratorVBR::make("gen", 0, 1, 2,
nomPower, nomPhPhVoltRMS, nomFreq, poleNum, nomFieldCurr,
Rs, Ll, Lmd, Lmd0, Lmq, Lmq0, Rfd, Llfd, Rkd, Llkd, Rkq1, Llkq1, Rkq2, Llkq2, H);
Real loadRes = 1037.8378;
Components::Base::Ptr r1 = Resistor::make("r1", 1, 0, loadRes);
Components::Base::Ptr r2 = Resistor::make("r2", 2, 0, loadRes);
Components::Base::Ptr r3 = Resistor::make("r3", 3, 0, loadRes);
Components::Base::Ptr r1 = Resistor::make("r1", 0, GND, loadRes);
Components::Base::Ptr r2 = Resistor::make("r2", 1, GND, loadRes);
Components::Base::Ptr r3 = Resistor::make("r3", 2, GND, loadRes);
Components::Base::List comps = { gen, r1, r2, r3 };
// Declare circuit components for resistance change
Real breakerRes = 0.001;
Components::Base::Ptr rBreaker1 = Resistor::make("rbreak1", 1, 0, breakerRes);
Components::Base::Ptr rBreaker2 = Resistor::make("rbreak2", 2, 0, breakerRes);
Components::Base::Ptr rBreaker3 = Resistor::make("rbreak3", 3, 0, breakerRes);
Components::Base::Ptr rBreaker1 = Resistor::make("rbreak1", 0, GND, breakerRes);
Components::Base::Ptr rBreaker2 = Resistor::make("rbreak2", 1, GND, breakerRes);
Components::Base::Ptr rBreaker3 = Resistor::make("rbreak3", 2, GND, breakerRes);
Components::Base::List compsBreakerOn = { rBreaker1, rBreaker2, rBreaker3, r1, r2, r3 };
Components::Base::List compsBreakerOn = {gen, rBreaker1, rBreaker2, rBreaker3, r1, r2, r3 };
// Set up simulation
Real tf, dt, t;
Real om = 2.0*M_PI*60.0;
tf = 0.3; dt = 0.000001; t = 0;
tf = 0.3; dt = 0.00001; t = 0;
Int downSampling = 50;
Simulation newSim("EMT_SynchronGenerator_VBR", comps, om, dt, tf, Logger::Level::INFO, SimulationType::EMT, downSampling);
newSim.setNumericalMethod(NumericalMethod::Trapezoidal_flux);
......
......@@ -9,8 +9,8 @@ Results_PLECS = csvread('../../../vsa/Results/SynGenDqEmt_ABCFault_Simulink/Volt
%omega_PLECS = csvread('../../../vsa/Results/SynGenDqEmt_ABCFault_Simulink/omega.csv');
%theta_PLECS = csvread('../../vsa/Results/SynGenVBREmt_ABCFault_PLECS/theta.csv');
%% read results from c++ simulation
VoltageVector = csvread('../../../vsa/Results/SynGenVbrEmt_ABCFault_DPsim/NewModel/data_vt.csv',1);
CurrentVector = csvread('../../../vsa/Results/SynGenVbrEmt_ABCFault_DPsim/NewModel/data_j.csv',1);
VoltageVector = csvread('../../../vsa/Results/SynGenVbrEmt_ABCFault_DPsim/NewModel/EMT_SynchronGenerator_VBR_LeftVector.csv',1);
CurrentVector = csvread('../../../vsa/Results/SynGenVbrEmt_ABCFault_DPsim/NewModel/EMT_SynchronGenerator_VBR_RightVector.csv',1);
%Log_SynGen = csvread('../../../vsa/Results/SynGenVbrEmt_ABCFault_DPsim/NewModel/SynGen_gen.csv',1);
%% Plot
figure(1)
......
......@@ -38,6 +38,7 @@ list(APPEND SOURCES
Components/EMT_VoltageSource_Ideal.cpp
Components/Exciter.cpp
Components/TurbineGovernor.cpp
Components/VoltageBehindReactanceEMTNew.cpp
)
find_package(Eigen3 REQUIRED)
......
/** Voltage behind reactance (EMT)
/** Voltage behind reactance (EMT)
*
* @author Markus Mirz <mmirz@eonerc.rwth-aachen.de>
* @copyright 2017, Institute for Automation of Complex Power Systems, EONERC
......
......@@ -24,47 +24,37 @@
using namespace DPsim;
VoltageBehindReactanceEMTNew::VoltageBehindReactanceEMTNew(String name, Int node1, Int node2, Int node3,
Components::EMT::VoltageBehindReactanceEMTNew::VoltageBehindReactanceEMTNew(String name, Int node1, Int node2, Int node3,
Real nomPower, Real nomVolt, Real nomFreq, Int poleNumber, Real nomFieldCur,
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, bool logActive)
: SynchGenBase(name, node1, node2, node3, nomPower, nomVolt, nomFreq, poleNumber, nomFieldCur,
Real inertia, Logger::Level logLevel)
: SynchronGeneratorBase(name, node1, node2, node3, nomPower, nomVolt, nomFreq, poleNumber, nomFieldCur,
Rs, Ll, Lmd, Lmd0, Lmq, Lmq0, Rfd, Llfd, Rkd, Llkd, Rkq1, Llkq1, Rkq2, Llkq2,
inertia, logActive) {
//mNumVirtualNodes = 3;
//mVirtualNodes = { 0, 0, 0 };
//mEa = IdealVoltageSourceEMT("mEa", 0, mVirtualNodes[0], mDVa);
//mEb = IdealVoltageSourceEMT("mEb", 0, mVirtualNodes[2], mDVb);
//mEc = IdealVoltageSourceEMT("mEc", 0, mVirtualNodes[4], mDVc);
inertia, logLevel)
{
}
VoltageBehindReactanceEMTNew::~VoltageBehindReactanceEMTNew() {
if (mLogActive) {
delete mLog;
}
Components::EMT::VoltageBehindReactanceEMTNew::~VoltageBehindReactanceEMTNew() {
}
void VoltageBehindReactanceEMTNew::AddExciter(Real Ta, Real Ka, Real Te, Real Ke, Real Tf, Real Kf, Real Tr, Real Lad, Real Rfd) {
void Components::EMT::VoltageBehindReactanceEMTNew::AddExciter(Real Ta, Real Ka, Real Te, Real Ke, Real Tf, Real Kf, Real Tr, Real Lad, Real Rfd) {
mExciter = Exciter(Ta, Ka, Te, Ke, Tf, Kf, Tr, Lad, Rfd);
mExciter.init(1, 1);
mExciter.initialize(1, 1);
WithExciter = true;
}
void VoltageBehindReactanceEMTNew::AddGovernor(Real Ta, Real Tb, Real Tc, Real Fa, Real Fb, Real Fc, Real K, Real Tsr, Real Tsm, Real Tm_init, Real PmRef) {
void Components::EMT::VoltageBehindReactanceEMTNew::AddGovernor(Real Ta, Real Tb, Real Tc, Real Fa, Real Fb, Real Fc, Real K, Real Tsr, Real Tsm, Real Tm_init, Real PmRef) {
mTurbineGovernor = TurbineGovernor(Ta, Tb, Tc, Fa, Fb, Fc, K, Tsr, Tsm);
mTurbineGovernor.init(PmRef, Tm_init);
mTurbineGovernor.initialize(PmRef, Tm_init);
WithTurbineGovernor = true;
}
void VoltageBehindReactanceEMTNew::init(Real om, Real dt,
void Components::EMT::VoltageBehindReactanceEMTNew::init(Real om, Real dt,
Real initActivePower, Real initReactivePower, Real initTerminalVolt, Real initVoltAngle, Real initFieldVoltage, Real initMechPower) {
mResistanceMat = Matrix::Zero(3, 3);
......@@ -173,7 +163,7 @@ void VoltageBehindReactanceEMTNew::init(Real om, Real dt,
}
void VoltageBehindReactanceEMTNew::step(SystemModel& system, Real time) {
void Components::EMT::VoltageBehindReactanceEMTNew::step(SystemModel& system, Real time) {
stepInPerUnit(system.getOmega(), system.getTimeStep(), time, system.getNumMethod());
......@@ -190,7 +180,8 @@ void VoltageBehindReactanceEMTNew::step(SystemModel& system, Real time) {
system.addRealToRightSideVector(mNode3, -mIc*mBase_i);
}
if (mLogActive) {
if (mLogLevel != Logger::Level::NONE) {
Matrix logValues(getRotorFluxes().rows() + getDqStatorCurrents().rows() + 3, 1);
logValues << getRotorFluxes(), getDqStatorCurrents(), getElectricalTorque(), getRotationalSpeed(), getRotorPosition();
mLog->LogDataLine(time, logValues);
......@@ -200,7 +191,7 @@ void VoltageBehindReactanceEMTNew::step(SystemModel& system, Real time) {
}
void VoltageBehindReactanceEMTNew::stepInPerUnit(Real om, Real dt, Real time, NumericalMethod numMethod) {
void Components::EMT::VoltageBehindReactanceEMTNew::stepInPerUnit(Real om, Real dt, Real time, NumericalMethod numMethod) {
mIabc <<
mIa,
......@@ -351,7 +342,7 @@ void VoltageBehindReactanceEMTNew::stepInPerUnit(Real om, Real dt, Real time, Nu
}
void VoltageBehindReactanceEMTNew::postStep(SystemModel& system) {
void Components::EMT::VoltageBehindReactanceEMTNew::postStep(SystemModel& system) {
if (mNode1 >= 0) {
mVa = system.getRealFromLeftSideVector(mNode1);
}
......@@ -372,7 +363,7 @@ void VoltageBehindReactanceEMTNew::postStep(SystemModel& system) {
}
}
void VoltageBehindReactanceEMTNew::CalculateLandpL() {
void Components::EMT::VoltageBehindReactanceEMTNew::CalculateLandpL() {
mDInductanceMat <<
mLl + mLa - mLb*cos(2 * mThetaMech), -mLa / 2 - mLb*cos(2 * mThetaMech - 2 * PI / 3), -mLa / 2 - mLb*cos(2 * mThetaMech + 2 * PI / 3),
-mLa / 2 - mLb*cos(2 * mThetaMech - 2 * PI / 3), mLl + mLa - mLb*cos(2 * mThetaMech - 4 * PI / 3), -mLa / 2 - mLb*cos(2 * mThetaMech),
......@@ -384,7 +375,7 @@ void VoltageBehindReactanceEMTNew::CalculateLandpL() {
pmDInductanceMat = pmDInductanceMat * 2 * mOmMech;
}
void VoltageBehindReactanceEMTNew::CalculateAuxiliarConstants(Real dt) {
void Components::EMT::VoltageBehindReactanceEMTNew::CalculateAuxiliarConstants(Real dt) {
b11 = (mRkq1 / mLlkq1)*(mDLmq / mLlkq1 - 1);
b12 = mRkq1*mDLmq / (mLlkq1*mLlkq2);
......@@ -447,7 +438,7 @@ void VoltageBehindReactanceEMTNew::CalculateAuxiliarConstants(Real dt) {
}
void VoltageBehindReactanceEMTNew::CalculateK() {
void Components::EMT::VoltageBehindReactanceEMTNew::CalculateK() {
c21_omega = -mOmMech*mDLmq / mLlkq1;
c22_omega = -mOmMech*mDLmq / mLlkq2;
......@@ -481,7 +472,7 @@ void VoltageBehindReactanceEMTNew::CalculateK() {
}
Matrix VoltageBehindReactanceEMTNew::parkTransform(Real theta, Real a, Real b, Real c) {
Matrix Components::EMT::VoltageBehindReactanceEMTNew::parkTransform(Real theta, Real a, Real b, Real c) {
Matrix dq0vector(3, 1);
......@@ -500,7 +491,7 @@ Matrix VoltageBehindReactanceEMTNew::parkTransform(Real theta, Real a, Real b, R
Matrix VoltageBehindReactanceEMTNew::inverseParkTransform(Real theta, Real q, Real d, Real zero) {
Matrix Components::EMT::VoltageBehindReactanceEMTNew::inverseParkTransform(Real theta, Real q, Real d, Real zero) {
Matrix abcVector(3, 1);
......
......@@ -22,246 +22,250 @@
#pragma once
#include "SynchGenBase.h"
#include "Base_SynchronGenerator.h"
#include "Exciter.h"
#include "TurbineGovernor.h"
#include "IdealvoltageSourceEMT.h"
#include "ResistorEMT.h"
#include "InductorEMT.h"
#include "EMT_VoltageSource_Ideal.h"
#include "EMT_Resistor.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 VoltageBehindReactanceEMTNew : public SynchGenBase {
protected:
/// Exciter Model
Exciter mExciter;
/// Determine if Exciter is activated
bool WithExciter = false;
/// Governor Model
TurbineGovernor mTurbineGovernor;
/// Determine if Turbine and Governor are activated
bool WithTurbineGovernor = false;
///// Voltage Source phase a
//IdealVoltageSourceEMT mEa;
///// Voltage Source phase b
//IdealVoltageSourceEMT mEb;
///// Voltage Source phase c
//IdealVoltageSourceEMT mEc;
/// d dynamic inductance
Real mDLmd;
/// q dynamic inductance
Real mDLmq;
/// Auxiliar inductance
Real mLa;
/// Auxiliar inductance
Real mLb;
/// d dynamic flux
Real mDPsid;
/// q dynamic flux
Real mDPsiq;
/// Dynamic d voltage
Real mDVq;
/// Dynamic q voltage
Real mDVd;
/// Dynamic voltage phase a
Real mDVa;
/// Dynamic voltage phase b
Real mDVb;
/// Dynamic voltage phase c
Real mDVc;
/// Interface voltage phase a
Real mVa;
/// Interface voltage phase b
Real mVb;
/// Interface voltage phase c
Real mVc;
/// Interface curent phase a
Real mIa;
/// Interface curent phase b
Real mIb;
/// Interface curent phase c
Real mIc;
/// Interface curent phase a last time step
Real mIa_hist;
/// Interface curent phase b last time step
Real mIb_hist;
/// Interface curent phase c last time step
Real mIc_hist;
///Load Resistance phase a
Real Ra;
///Load Resistance phase b
Real Rb;
///Load Resistance phase c
Real Rc;
/// Magnetizing flux linkage in q axis
Real mPsimq;
/// Magnetizing flux linkage in d axis
Real mPsimd;
/// Rotor flux vector
Matrix mRotorFlux = Matrix::Zero(4, 1);
/// Dq stator current vector
Matrix mDqStatorCurrents = Matrix::Zero(2, 1);
/// Dq stator current vector - from previous time step
Matrix mDqStatorCurrents_hist = Matrix::Zero(2, 1);
// ### Useful Matrices ###
/// inductance matrix
Matrix mDInductanceMat = Matrix::Zero(3, 3);
/// Derivative of inductance matrix
Matrix pmDInductanceMat = Matrix::Zero(3, 3);
/// Load Resistance
Matrix R_load = Matrix::Zero(3, 3);
/// Phase currents in pu
Matrix mIabc = Matrix::Zero(3, 1);
/// Subtransient voltage in pu
Matrix mDVabc = Matrix::Zero(3, 1);
Matrix mDVabc_hist = Matrix::Zero(3, 1);
/// Matrix paremeters for integration of rotor flux linkages - A
Matrix A_flux = Matrix::Zero(4, 4);
/// Variables for integration of rotor flux linkages - B
Matrix B_flux = Matrix::Zero(4, 2);
/// Variables for integration of rotor flux linkages - C
Matrix C_flux = Matrix::Zero(4, 1);
/// Auxiliar variables
Real c21_omega;
Real c22_omega;
Real c13_omega;
Real c14_omega;
Matrix K1a = Matrix::Zero(2, 2);
Matrix K1b = Matrix::Zero(2, 1);
Matrix K1 = Matrix::Zero(2, 1);
Matrix K2a = Matrix::Zero(2, 2);
Matrix K2b = Matrix::Zero(2, 1);
Matrix K2 = Matrix::Zero(2, 1);
Matrix H_qdr;
Matrix K = Matrix::Zero(3, 3);
Matrix mPsikq1kq2 = Matrix::Zero(2, 1);
Matrix mPsifdkd = Matrix::Zero(2, 1);
/// Auxiliar constants
Real c11;
Real c12;
Real c23;
Real c24;
Real c15;
Real c25;
Real c26;
Real b11;
Real b12;
Real b13;
Real b21;
Real b22;
Real b23;
Real b31;
Real b32;
Real b33;
Real b41;
Real b42;
Real b43;
Matrix Ea = Matrix::Zero(2, 2);
Matrix E1b = Matrix::Zero(2, 1);
Matrix E1 = Matrix::Zero(2, 2);
Matrix Fa = Matrix::Zero(2, 2);
Matrix F1b = Matrix::Zero(2, 1);
Matrix F1 = Matrix::Zero(2, 2);
Matrix E2b = Matrix::Zero(2, 2);
Matrix E2 = Matrix::Zero(2, 2);
Matrix F2b = Matrix::Zero(2, 2);
Matrix F2 = Matrix::Zero(2, 2);
Matrix F3b = Matrix::Zero(2, 1);
Matrix F3 = Matrix::Zero(2, 2);
Matrix C26 = Matrix::Zero(2, 1);
public:
VoltageBehindReactanceEMTNew() { };
~VoltageBehindReactanceEMTNew();
/// 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.
VoltageBehindReactanceEMTNew(String name, Int node1, Int node2, Int node3,
Real nomPower, Real nomVolt, Real nomFreq, Int poleNumber, Real nomFieldCur,
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, bool logActive = false);
/// Function to initialize Exciter
void AddExciter(Real Ta, Real Ka, Real Te, Real Ke, Real Tf, Real Kf, Real Tr, Real Lad, Real Rfd);
/// Function to initialize Governor and Turbine
void AddGovernor(Real Ta, Real Tb, Real Tc, Real Fa, Real Fb, Real Fc, Real K, Real Tsr, Real Tsm, Real Tm_init, Real PmRef);
/// 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, Real initFieldVoltage, Real initMechPower);
/// 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 time);
/// 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 time, NumericalMethod numMethod);
/// Retrieves calculated voltage from simulation for next step
void postStep(SystemModel& system);
/// Park transform as described in Krause
Matrix parkTransform(Real theta, Real a, Real b, Real c);
/// Inverse Park transform as described in Krause
Matrix inverseParkTransform(Real theta, Real q, Real d, Real zero);
/// Calculate inductance Matrix L and its derivative
void CalculateLandpL();
void CalculateAuxiliarConstants(Real dt);
void CalculateK();
Matrix& getRotorFluxes() { return mRotorFlux; }
Matrix& getDqStatorCurrents() { return mDqStatorCurrents; }
Real getElectricalTorque() { return mElecTorque*mBase_T; }
Real getRotationalSpeed() { return mOmMech*mBase_OmMech; }
Real getRotorPosition() { return mThetaMech; }
void init(Real om, Real dt) { }
void applySystemMatrixStamp(SystemModel& system) { }
void applyRightSideVectorStamp(SystemModel& system) { }
};
namespace DPsim {
namespace Components {
namespace EMT {
/// 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 VoltageBehindReactanceEMTNew : public SynchronGeneratorBase, public SharedFactory<VoltageBehindReactanceEMTNew> {
protected:
/// Exciter Model
Exciter mExciter;
/// Determine if Exciter is activated
bool WithExciter = false;
/// Governor Model
TurbineGovernor mTurbineGovernor;
/// Determine if Turbine and Governor are activated
bool WithTurbineGovernor = false;
///// Voltage Source phase a
//IdealVoltageSourceEMT mEa;
///// Voltage Source phase b
//IdealVoltageSourceEMT mEb;
///// Voltage Source phase c
//IdealVoltageSourceEMT mEc;
/// d dynamic inductance
Real mDLmd;
/// q dynamic inductance
Real mDLmq;
/// Auxiliar inductance
Real mLa;
/// Auxiliar inductance
Real mLb;
/// d dynamic flux
Real mDPsid;
/// q dynamic flux
Real mDPsiq;
/// Dynamic d voltage
Real mDVq;
/// Dynamic q voltage
Real mDVd;
/// Dynamic voltage phase a
Real mDVa;
/// Dynamic voltage phase b
Real mDVb;
/// Dynamic voltage phase c
Real mDVc;
/// Interface voltage phase a
Real mVa;
/// Interface voltage phase b
Real mVb;
/// Interface voltage phase c
Real mVc;
/// Interface curent phase a
Real mIa;
/// Interface curent phase b
Real mIb;
/// Interface curent phase c
Real mIc;
/// Interface curent phase a last time step
Real mIa_hist;
/// Interface curent phase b last time step
Real mIb_hist;
/// Interface curent phase c last time step
Real mIc_hist;
///Load Resistance phase a
Real Ra;
///Load Resistance phase b
Real Rb;
///Load Resistance phase c
Real Rc;
/// Magnetizing flux linkage in q axis
Real mPsimq;
/// Magnetizing flux linkage in d axis
Real mPsimd;
/// Rotor flux vector
Matrix mRotorFlux = Matrix::Zero(4, 1);
/// Dq stator current vector
Matrix mDqStatorCurrents = Matrix::Zero(2, 1);
/// Dq stator current vector - from previous time step
Matrix mDqStatorCurrents_hist = Matrix::Zero(2, 1);
// ### Useful Matrices ###
/// inductance matrix
Matrix mDInductanceMat = Matrix::Zero(3, 3);
/// Derivative of inductance matrix
Matrix pmDInductanceMat = Matrix::Zero(3, 3);
/// Load Resistance
Matrix R_load = Matrix::Zero(3, 3);
/// Phase currents in pu