Commit 66a41836 authored by Viviane Sapucaia's avatar Viviane Sapucaia
Browse files

Fixed error in Exciter model

Added exciter model to detailed VBR model
parent 49184cc1
......@@ -81,7 +81,7 @@ int main() {
circElements.push_back(r3);
// Declare circuit components for resistance change
Real breakerRes = 0.001;
Real breakerRes = 103.78378;
ElementPtr rBreaker1 = std::make_shared<ResistorEMT>("rbreak1", 1, 0, breakerRes);
ElementPtr rBreaker2 = std::make_shared<ResistorEMT>("rbreak2", 2, 0, breakerRes);
ElementPtr rBreaker3 = std::make_shared<ResistorEMT>("rbreak3", 3, 0, breakerRes);
......@@ -90,14 +90,14 @@ int main() {
circElementsBreakerOn.push_back(rBreaker1);
circElementsBreakerOn.push_back(rBreaker2);
circElementsBreakerOn.push_back(rBreaker3);
circElementsBreakerOn.push_back(r1);
circElementsBreakerOn.push_back(r2);
circElementsBreakerOn.push_back(r3);
//circElementsBreakerOn.push_back(r1);
//circElementsBreakerOn.push_back(r2);
//circElementsBreakerOn.push_back(r3);
// Set up simulation
Real tf, dt, t;
Real om = 2.0*M_PI*60.0;
tf = 0.3; dt = 0.0001; t = 0;
tf = 10; dt = 0.0001; t = 0;
Int downSampling = 1;
Simulation newSim(circElements, om, dt, tf, log, SimulationType::EMT, downSampling);
newSim.setNumericalMethod(NumericalMethod::Trapezoidal_flux);
......@@ -113,7 +113,7 @@ int main() {
Real mechPower = 5.5558e5;
shared_ptr<VoltageBehindReactanceEMT> genPtr = std::dynamic_pointer_cast<VoltageBehindReactanceEMT>(gen);
genPtr->init(om, dt, initActivePower, initReactivePower, initTerminalVolt, initVoltAngle, fieldVoltage, mechPower);
//genPtr->AddExciter(Ta, Ka, Te, Ke, Tf, Kf, Tr, Lmd, Rfd);
genPtr->AddExciter(Ta, Ka, Te, Ke, Tf, Kf, Tr, Lmd, Rfd);
// Calculate initial values for circuit at generator connection point
Real initApparentPower = sqrt(pow(initActivePower, 2) + pow(initReactivePower, 2));
......@@ -129,8 +129,8 @@ int main() {
Real lastLogTime = 0;
Real logTimeStep = 0.0001;
newSim.setSwitchTime(0.1, 1);
newSim.setSwitchTime(0.2, 0);
newSim.setSwitchTime(2, 1);
newSim.setSwitchTime(6, 0);
// Main Simulation Loop
while (newSim.getTime() < tf) {
......
......@@ -3,9 +3,9 @@
%% read PLECS results
Results_PLECS = csvread('../../../vsa/Results/SynGenVBREmt_ABCFault_PLECS/Voltages_and_currents.csv');
%Te_PLECS = csvread('../../vsa/Results/SynGenVBREmt_ABCFault_PLECS/electrical_torque.csv');
omega_PLECS = csvread('../../../vsa/Results/SynGenVBREmt_ABCFault_PLECS/omega.csv');
Results_PLECS = csvread('../../../vsa/Results/SimpSynGenDqEmtExciter_LoadChange_Simulink/Voltages_and_currents.csv');
Te_PLECS = csvread('../../../vsa/Results/SimpSynGenDqEmtExciter_LoadChange_Simulink/electrical_torque.csv');
omega_PLECS = csvread('../../../vsa/Results/SimpSynGenDqEmtExciter_LoadChange_Simulink/omega.csv');
%theta_PLECS = csvread('../../vsa/Results/SynGenVBREmt_ABCFault_PLECS/theta.csv');
%% read results from c++ simulation
VoltageVector = csvread('../../../vsa/Results/Testing/data_vt.csv');
......@@ -16,31 +16,31 @@ figure(1)
hold off
plot(VoltageVector(:,1),VoltageVector(:,2));
hold on
plot(VoltageVector(:,1),VoltageVector(:,3));
plot(VoltageVector(:,1),VoltageVector(:,4));
%plot(VoltageVector(:,1),VoltageVector(:,3));
%plot(VoltageVector(:,1),VoltageVector(:,4));
plot(Results_PLECS(:,1),Results_PLECS(:,2),'--');
plot(Results_PLECS(:,1),Results_PLECS(:,3),'--');
plot(Results_PLECS(:,1),Results_PLECS(:,4),'--');
%plot(Results_PLECS(:,1),Results_PLECS(:,3),'--');
%plot(Results_PLECS(:,1),Results_PLECS(:,4),'--');
title('Phase Voltages');
legend('va DPSim','vb DPSim', 'vc DPSim','va PLECS','vb PLECS','vc PLECS');
%legend('va DPSim','vb DPSim', 'vc DPSim','va Simulink','vb Simulink','vc Simulink');
%legend('va DPSim','vb DPSim', 'vc DPSim','va PLECS','vb PLECS','vc PLECS');
legend('va DPSim','va Simulink');
figure(2)
hold off
plot(CurrentVector(:,1),-CurrentVector(:,2));
plot(CurrentVector(:,1),CurrentVector(:,2));
hold on
plot(CurrentVector(:,1),-CurrentVector(:,3));
plot(CurrentVector(:,1),-CurrentVector(:,4));
%plot(CurrentVector(:,1),CurrentVector(:,3));
%plot(CurrentVector(:,1),CurrentVector(:,4));
plot(Results_PLECS(:,1),Results_PLECS(:,5),'--');
plot(Results_PLECS(:,1),Results_PLECS(:,6),'--');
plot(Results_PLECS(:,1),Results_PLECS(:,7),'--');
%plot(Results_PLECS(:,1),Results_PLECS(:,6),'--');
%plot(Results_PLECS(:,1),Results_PLECS(:,7),'--');
title('Phase Currents');
legend('ia DPSim','ib DPSim','ic DPSim','ia PLECS','ib PLECS','ic PLECS');
%legend('ia DPSim','ib DPSim','ic DPSim','ia Simulink','ib Simulink','ic Simulink');
%legend('ia DPSim','ib DPSim','ic DPSim','ia PLECS','ib PLECS','ic PLECS');
legend('ia DPSim','ia Simulink');
% figure(3)
% hold off
......@@ -78,14 +78,14 @@ plot(Results_PLECS(:,1),omega_PLECS);
title('Rotor speed');
legend('\omega DPSim','\omega PLECS');
%
% figure(4)
% hold off
% plot(Log_SynGen(:,1),Log_SynGen(:,7));
% hold on
% plot(Results_PLECS(:,1),Te_PLECS);
%
% title('Electrical Torque');
% legend('Te DPSim','Te PLECS');
figure(4)
hold off
plot(Log_SynGen(:,1),Log_SynGen(:,7));
hold on
plot(Results_PLECS(:,1),-Te_PLECS);
title('Electrical Torque');
legend('Te DPSim','Te PLECS');
%
% figure(5)
% hold off
......
......@@ -37,30 +37,25 @@ Exciter::Exciter(Real Ta, Real Ka, Real Te, Real Ke, Real Tf, Real Kf, Real Tr,
mTr = Tr;
mLad = Lad;
mRfd = Rfd;
mVf = Vf_init*(Lad / Rfd);
mVf = Vf_init*(mLad / mRfd);
}
Exciter::~Exciter() {
}
Real Exciter::step(Real mVd, Real mVq, Real Vref, Real dt) {
mVh = sqrt(pow(mVd, 2.) + pow(mVq, 2.));
// Voltage Transducer equation
mVm = Euler(mVm, -1, 1, dt / mTr, mVh);
mVm = Trapezoidal(mVm, -1, 1, dt / mTr, mVh);
// Stabilizing feedback equation
Real dUf = (mVr - mVse - mKe*mVf) / mTe;
mVis = Euler(mVis, -1, mKf, dt / mTf, dUf);
mVfl = mVfl + dt*mVis / mTf;
mVis = (mKf / mTf)*mVf - mVfl;
// Amplifier equation
mVr = Euler(mVr, -1, mKa, dt / mTa, Vref - mVm - mVis);
mVr = Trapezoidal(mVr, -1, mKa, dt / mTa, Vref - mVm - mVis);
// Exciter
mVse = 0.0039*exp(mVf*1.555);
mVf = Euler(mVf, -mKe, 1, dt / mTe, mVr - mVse);
mVf = Trapezoidal(mVf, -mKe, 1, dt / mTe, mVr - mVse);
return (mRfd / mLad)*mVf;
......
......@@ -48,21 +48,22 @@ namespace DPsim {
Real mLad;
Real mRfd;
// Reference voltage
/// Reference voltage
Real mVref = 0;
// Output of voltage transducer
/// Output of voltage transducer
Real mVm = 0;
// Input of voltage transducer
/// Input of voltage transducer
Real mVh = 0;
// Output of stablizing feedback
/// Output of stablizing feedback
Real mVis = 0;
// Output of se function
/// Output of se function
Real mVse = 0;
// Regulator output
/// Regulator output
Real mVr = 0;
// Exciter output
/// Exciter output
Real mVf = 0;
Real mVfd;
/// Auxiliar variable
Real mVfl = 0;
bool mLogActive;
Logger* mLog;
......@@ -70,7 +71,7 @@ namespace DPsim {
public:
Exciter() { };
~Exciter();
~Exciter() {};
/// Initializes exciter parameters
Exciter(Real Ta, Real Ka, Real Te, Real Ke, Real Tf, Real Kf, Real Tr, Real Lad, Real Rfd, Real Vf_init);
......@@ -79,7 +80,7 @@ namespace DPsim {
Real step(Real mVd, Real mVq, Real Vref, Real dt);
Real UpdateFieldVoltage() { return mVfd; }
};
}
......@@ -125,6 +125,8 @@ void SimplifiedVBR::init(Real om, Real dt,
void SimplifiedVBR::step(SystemModel& system, Real time) {
R_load = system.getCurrentSystemMatrix().inverse() / mBase_Z;
stepInPerUnit(system.getOmega(), system.getTimeStep(), time, system.getNumMethod());
// Update current source accordingly
......@@ -150,7 +152,8 @@ void SimplifiedVBR::step(SystemModel& system, Real time) {
void SimplifiedVBR::stepInPerUnit(Real om, Real dt, Real time, NumericalMethod numMethod) {
// Calculate mechanical variables with euler
mElecTorque = (mPsimd*mIq - mPsimq*mId);
//mElecTorque = (mPsimd*mIq - mPsimq*mId);
mElecTorque = mVd*mId + mVq*mIq + mRs*(mId*mId + mIq*mIq);
mOmMech = mOmMech + dt * (1. / (2. * mH) * (mElecTorque - mMechTorque));
mThetaMech = mThetaMech + dt * (mOmMech* mBase_OmMech);
......@@ -172,9 +175,10 @@ void SimplifiedVBR::stepInPerUnit(Real om, Real dt, Real time, NumericalMethod n
mIb = inverseParkTransform(mThetaMech, mIq, mId, 0)(1);
mIc = inverseParkTransform(mThetaMech, mIq, mId, 0)(2);
mVq = -mRs*mIq - mOmMech*mDLd*mId + mDVq;
mVd = -mRs*mId + mOmMech*mDLq*mIq + mDVd;
if (WithExciter == true) {
mVq = -mRs*mIq - mOmMech*mDLd*mId + mDVq;
mVd = -mRs*mId + mOmMech*mDLq*mIq + mDVd;
mVfd = mExciter.step(mVd, mVq, 1, dt);
}
......@@ -215,44 +219,6 @@ void SimplifiedVBR::stepInPerUnit(Real om, Real dt, Real time, NumericalMethod n
mDVd;
// Load resistance
if (time < 0.1)
{
Ra = 1037.8378 / mBase_Z;
Rb = 1037.8378 / mBase_Z;
Rc = 1037.8378 / mBase_Z;
R_load <<
Ra, 0, 0,
0, Rb, 0,
0, 0, Rc;
}
else if (time > 0.2)
{
if (signbit(mIa) != signbit(mIa_hist))
Ra = 1037.8378 / mBase_Z;
if (signbit(mIb) != signbit(mIb_hist))
Rb = 1037.8378 / mBase_Z;
if (signbit(mIc) != signbit(mIc_hist))
Rc = 1037.8378 / mBase_Z;
R_load <<
Ra, 0, 0,
0, Rb, 0,
0, 0, Rc;
}
else
{
Ra = 0.001 / mBase_Z;
Rb = 0.001 / mBase_Z;
Rc = 0.001 / mBase_Z;
R_load <<
Ra, 0, 0,
0, Rb, 0,
0, 0, Rc;
}
}
......
......@@ -43,6 +43,11 @@ VoltageBehindReactanceEMT::~VoltageBehindReactanceEMT() {
}
}
void VoltageBehindReactanceEMT::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, mVfd);
WithExciter = true;
}
void VoltageBehindReactanceEMT::init(Real om, Real dt,
......@@ -153,9 +158,6 @@ void VoltageBehindReactanceEMT::init(Real om, Real dt,
void VoltageBehindReactanceEMT::step(SystemModel& system, Real time) {
stepInPerUnit(system.getOmega(), system.getTimeStep(), time, system.getNumMethod());
R_load = system.getCurrentSystemMatrix().inverse() / mBase_Z;
......@@ -221,6 +223,15 @@ void VoltageBehindReactanceEMT::stepInPerUnit(Real om, Real dt, Real time, Numer
mId = parkTransform(mThetaMech, mIa, mIb, mIc)(1);
mI0 = parkTransform(mThetaMech, mIa, mIb, mIc)(2);
if (WithExciter == true) {
// dq-transform of interface voltage
mVd = parkTransform(mThetaMech, mVa / mBase_v, mVb / mBase_v, mVc / mBase_v)(0);
mVq = parkTransform(mThetaMech, mVa / mBase_v, mVb / mBase_v, mVc / mBase_v)(1);
mV0 = parkTransform(mThetaMech, mVa / mBase_v, mVb / mBase_v, mVc / mBase_v)(2);
mVfd = mExciter.step(mVd, mVq, 1, dt);
}
// Calculate rotor flux likanges
if (mNumDampingWindings == 2)
{
......@@ -326,7 +337,24 @@ void VoltageBehindReactanceEMT::stepInPerUnit(Real om, Real dt, Real time, Numer
void VoltageBehindReactanceEMT::postStep(SystemModel& system) {
if (mNode1 >= 0) {
mVa = system.getRealFromLeftSideVector(mNode1);
}
else {
mVa = 0;
}
if (mNode2 >= 0) {
mVb = system.getRealFromLeftSideVector(mNode2);
}
else {
mVb = 0;
}
if (mNode3 >= 0) {
mVc = system.getRealFromLeftSideVector(mNode3);
}
else {
mVc = 0;
}
}
void VoltageBehindReactanceEMT::CalculateLandpL() {
......@@ -346,14 +374,15 @@ Matrix VoltageBehindReactanceEMT::parkTransform(Real theta, Real a, Real b, Real
Matrix dq0vector(3, 1);
Real q, d;
Real q, d, zero;
q = 2. / 3. * cos(theta)*a + 2. / 3. * cos(theta - 2. * M_PI / 3.)*b + 2. / 3. * cos(theta + 2. * M_PI / 3.)*c;
d = 2. / 3. * sin(theta)*a + 2. / 3. * sin(theta - 2. * M_PI / 3.)*b + 2. / 3. * sin(theta + 2. * M_PI / 3.)*c;
zero = 1. / 3.*a + 1. / 3.*b + 1. / 3.*c;
dq0vector << q,
d,
0;
zero;
return dq0vector;
}
......
......@@ -24,6 +24,7 @@
#pragma once
#include "SynchGenBase.h"
#include "Exciter.h"
namespace DPsim {
......@@ -37,6 +38,10 @@ namespace DPsim {
class VoltageBehindReactanceEMT : public SynchGenBase {
protected:
/// Exciter Model
Exciter mExciter;
bool WithExciter = false;
/// d dynamic inductance
Real mDLmd;
/// q dynamic inductance
......@@ -142,6 +147,8 @@ namespace DPsim {
Real Rkq1, Real Llkq1, Real Rkq2, Real Llkq2,
Real inertia, bool logActive = false);
void AddExciter(Real Ta, Real Ka, Real Te, Real Ke, Real Tf, Real Kf, Real Tr, Real Lad, Real Rfd);
/// 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,
......
......@@ -202,15 +202,16 @@ Int Simulation::stepGeneratorTest(Logger& leftSideVectorLog, Logger& rightSideVe
if (mCurrentSwitchTimeIndex < mSwitchEventVector.size()) {
if (mTime >= mSwitchEventVector[mCurrentSwitchTimeIndex].switchTime) {
if (mCurrentSwitchTimeIndex == 1) {
clearFault(1, 2, 3);
mCurrentSwitchTimeIndex++;
}
else {
switchSystemMatrix(mSwitchEventVector[mCurrentSwitchTimeIndex].systemIndex);
mElements = mElementsVector[mSwitchEventVector[mCurrentSwitchTimeIndex++].systemIndex];
}
switchSystemMatrix(mSwitchEventVector[mCurrentSwitchTimeIndex].systemIndex);
mElements = mElementsVector[mSwitchEventVector[mCurrentSwitchTimeIndex++].systemIndex];
//if (mCurrentSwitchTimeIndex == 1) {
// clearFault(1, 2, 3);
// mCurrentSwitchTimeIndex++;
//}
//else {
// switchSystemMatrix(mSwitchEventVector[mCurrentSwitchTimeIndex].systemIndex);
// mElements = mElementsVector[mSwitchEventVector[mCurrentSwitchTimeIndex++].systemIndex];
//}
//mCurrentSwitchTimeIndex++;
mLogger->Log(LogLevel::INFO) << "Switched to system " << mCurrentSwitchTimeIndex << " at " << mTime << std::endl;
mLogger->Log(LogLevel::INFO) << "New matrix:" << std::endl << mSystemModel.getCurrentSystemMatrix() << std::endl;
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment