Commit 53fcbc1d authored by Jan Dinkelbach's avatar Jan Dinkelbach Committed by Markus Mirz
Browse files

refactor task definition of trafo in sp1ph and add validation notebook for trafo

parent 58bfd91a
......@@ -41,6 +41,7 @@ set(CIRCUIT_SOURCES
# SP examples
Circuits/SP_Circuits.cpp
Circuits/SP_PiLine.cpp
Circuits/SP_Trafo.cpp
# Combined EMT/DP/SP examples
Circuits/EMT_DP_SP_Slack.cpp
......
/* Copyright 2017-2020 Institute for Automation of Complex Power Systems,
* EONERC, RWTH Aachen University
* DPsim
*
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* any later version.
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
* You should have received a copy of the GNU General Public License
* along with this program. If not, see <http://www.gnu.org/licenses/>.
*********************************************************************************/
#include <DPsim.h>
using namespace DPsim;
using namespace CPS;
void simTrafoElementsSP1ph() {
Real timeStep = 0.00005;
Real finalTime = 1;
String simName = "SP_Trafo_Elements";
Logger::setLogDir("logs/"+simName);
Real voltageHVSide = 100000;
Real voltageMVSide = 10000;
Real trafoResistance = 1;
Real trafoInductance = 0.1;
Real loadResistanceHVSide = 10000;
Real ratio = voltageHVSide/voltageMVSide;
Real snubberResistanceHVSide = ratio*ratio*voltageMVSide*1e6;
// Nodes
auto n1 = SimNode<Complex>::make("n1");
auto n2 = SimNode<Complex>::make("n2");
auto vn1 = SimNode<Complex>::make("vn1");
// Components
auto vs = SP::Ph1::VoltageSource::make("v_1");
auto trafoRes = SP::Ph1::Resistor::make("trafo_res");
auto trafoSnubberRes = SP::Ph1::Resistor::make("trafo_snub_res");
auto trafoInd = SP::Ph1::Inductor::make("trafo_ind");
auto loadRes = SP::Ph1::Resistor::make("r_1");
// Topology
vs->connect({ SimNode<Complex>::GND, n1 });
trafoRes->connect({ n1, vn1 });
trafoInd->connect({ vn1, n2 });
trafoSnubberRes->connect({ n2, SimNode<Complex>::GND });
loadRes->connect({ n2, SimNode<Complex>::GND });
// Parameters
vs->setParameters(CPS::Math::polar(voltageHVSide, 0));
trafoRes->setParameters(trafoResistance);
trafoInd->setParameters(trafoInductance);
trafoSnubberRes->setParameters(snubberResistanceHVSide);
loadRes->setParameters(loadResistanceHVSide);
// Define system topology
SystemTopology sys(50, SystemNodeList{n1, n2, vn1 }, SystemComponentList{vs, trafoRes, trafoInd, trafoSnubberRes, loadRes});
// Logging
auto logger = DataLogger::make(simName);
logger->addAttribute("v1", n1->attribute("v"));
logger->addAttribute("itrafo", trafoInd->attribute("i_intf"));
Simulation sim(simName);
sim.setSystem(sys);
sim.setTimeStep(timeStep);
sim.setFinalTime(finalTime);
sim.setDomain(Domain::SP);
sim.addLogger(logger);
sim.run();
}
void simTrafoSP1ph() {
Real timeStep = 0.00005;
Real finalTime = 1;
String simName = "SP_Trafo_Component";
Logger::setLogDir("logs/"+simName);
Real voltageHVSide = 100000;
Real voltageMVSide = 10000;
Real trafoResistance = 1;
Real trafoInductance = 0.1;
Real loadResistanceHVSide = 10000;
Real ratio = voltageHVSide/voltageMVSide;
Real loadResistanceMVSide = loadResistanceHVSide/(ratio*ratio);
// Nodes
auto n1 = SimNode<Complex>::make("n1");
auto n2 = SimNode<Complex>::make("n2");
// Components
auto vs = SP::Ph1::VoltageSource::make("v_1", Logger::Level::debug);
auto trafo = std::make_shared<SP::Ph1::Transformer>("trafo", "trafo", Logger::Level::debug, true);
auto loadRes = SP::Ph1::Resistor::make("r_1", Logger::Level::debug);
// Topology
vs->connect({ SimNode<Complex>::GND, n1 });
trafo->connect({ n1, n2 });
loadRes->connect({ n2, SimNode<Complex>::GND });
// Parameters
vs->setParameters(CPS::Math::polar(voltageHVSide, 0));
trafo->setParameters(voltageHVSide, voltageMVSide, 50e6, ratio, 0, trafoResistance, trafoInductance, 2.0*PI*50.0);
loadRes->setParameters(loadResistanceMVSide);
// Define system topology
SystemTopology sys(50, SystemNodeList{n1, n2 }, SystemComponentList{vs, trafo, loadRes});
// Logging
auto logger = DataLogger::make(simName);
logger->addAttribute("v1", n1->attribute("v"));
logger->addAttribute("itrafo", trafo->attribute("i_intf"));
Simulation sim(simName);
sim.setSystem(sys);
sim.setTimeStep(timeStep);
sim.setFinalTime(finalTime);
sim.setDomain(Domain::SP);
sim.addLogger(logger);
sim.run();
}
int main(int argc, char* argv[]) {
simTrafoElementsSP1ph();
simTrafoSP1ph();
return 0;
}
%% Cell type:markdown id: tags:
# Trafo Tests
%% Cell type:code id: tags:
``` python
import villas.dataprocessing.readtools as rt
from villas.dataprocessing.timeseries import TimeSeries as ts
import matplotlib.pyplot as plt
import numpy as np
# %matplotlib widget
epsilon = 1e-12
```
%% Cell type:markdown id: tags:
## SP Trafo with elements
%% Cell type:code id: tags:
``` python
%%bash
TOP=${TOP:-$(git rev-parse --show-toplevel)}
PATH=${TOP}/build/Examples/Cxx
SP_Trafo
```
%% Cell type:code id: tags:
``` python
work_dir = 'logs/SP_Trafo_Elements/'
log_name = 'SP_Trafo_Elements'
print(work_dir + log_name + '.csv')
trafo_elements = rt.read_timeseries_dpsim(work_dir + log_name + '.csv')
trafo_elements_sp_shifted = ts.frequency_shift_list(trafo_elements, 50)
```
%% Cell type:code id: tags:
``` python
plt.figure()
plt.plot(trafo_elements_sp_shifted['v1_shift'].time, trafo_elements_sp_shifted['v1_shift'].values, label='v1_shift')
plt.legend()
```
%% Cell type:code id: tags:
``` python
plt.figure()
plt.plot(trafo_elements_sp_shifted['itrafo_shift'].time, trafo_elements_sp_shifted['itrafo_shift'].values, label='itrafo_shift')
plt.legend()
```
%% Cell type:markdown id: tags:
## SP Trafo composite model
%% Cell type:code id: tags:
``` python
work_dir = 'logs/SP_Trafo_Component/'
log_name = 'SP_Trafo_Component'
print(work_dir + log_name + '.csv')
trafo_component = rt.read_timeseries_dpsim(work_dir + log_name + '.csv')
trafo_component_sp_shifted = ts.frequency_shift_list(trafo_component, 50)
```
%% Cell type:code id: tags:
``` python
plt.figure()
plt.plot(trafo_component_sp_shifted['v1_shift'].time, trafo_component_sp_shifted['v1_shift'].values, label='v1_shift')
plt.legend()
```
%% Cell type:code id: tags:
``` python
plt.figure()
plt.plot(trafo_component_sp_shifted['itrafo_shift'].time, trafo_component_sp_shifted['itrafo_shift'].values, label='itrafo_shift')
plt.legend()
```
%% Cell type:markdown id: tags:
## Error for SP Trafo
%% Cell type:code id: tags:
``` python
plt.figure()
for name in ['v1_shift', 'itrafo_shift']:
plt.plot(trafo_elements_sp_shifted[name].time, trafo_elements_sp_shifted[name].values - trafo_component_sp_shifted[name].values, label=name+'_error')
plt.legend()
```
%% Cell type:markdown id: tags:
## Assertion for SP Trafo
%% Cell type:code id: tags:
``` python
errors_sp_shifted = []
for name in ['v1_shift', 'itrafo_shift']:
errors_sp_shifted.append(np.absolute(trafo_elements_sp_shifted[name].values - trafo_component_sp_shifted[name].values).max())
print(name + ': ' + str(errors_sp_shifted[-1]))
assert np.max(errors_sp_shifted) < epsilon
```
%% Cell type:code id: tags:
``` python
```
/*
* @copyright 2017, Institute for Automation of Complex Power Systems, EONERC
/* Copyright 2017-2020 Institute for Automation of Complex Power Systems,
* EONERC, RWTH Aachen University
*
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program. If not, see <http://www.gnu.org/licenses/>.
* This Source Code Form is subject to the terms of the Mozilla Public
* License, v. 2.0. If a copy of the MPL was not distributed with this
* file, You can obtain one at https://mozilla.org/MPL/2.0/.
*********************************************************************************/
#pragma once
#include <cps/SimPowerComp.h>
#include <cps/Solver/PFSolverInterfaceBranch.h>
#include <cps/Solver/MNAInterface.h>
......@@ -43,6 +35,9 @@ namespace Ph1 {
///
std::shared_ptr<SP::Ph1::Resistor> mSubResistor;
/// Snubber resistance added on the low voltage side
Real mSnubberResistance;
/// Rated Apparent Power [VA]
Real mRatedPower = 0;
/// Transformer ratio magnitude
......@@ -115,10 +110,12 @@ namespace Ph1 {
/// nodal reactive power injection
Real mReactivePowerInjection;
/// Boolean for considering resistive losses with sub resistor
Bool mWithResistiveLosses;
public:
/// Defines UID, name and logging level
Transformer(String uid, String name,
Logger::Level logLevel = Logger::Level::off);
Logger::Level logLevel = Logger::Level::off, Bool withResistiveLosses = false);
/// Defines name and logging level
Transformer(String name, Logger::Level logLevel = Logger::Level::off)
: Transformer(name, name, logLevel) { }
......@@ -156,18 +153,17 @@ namespace Ph1 {
void mnaUpdateCurrent(const Matrix& leftVector) override;
/// Updates internal voltage variable of the component
void mnaUpdateVoltage(const Matrix& leftVector) override;
/// MNA post step operations
void mnaPostStep(Real time, Int timeStepCount, Attribute<Matrix>::Ptr &leftVector);
/// Add MNA post step dependencies
void mnaAddPostStepDependencies(AttributeBase::List &prevStepDependencies, AttributeBase::List &attributeDependencies, AttributeBase::List &modifiedAttributes, Attribute<Matrix>::Ptr &leftVector);
class MnaPostStep : public Task {
public:
MnaPostStep(Transformer& transformer, Attribute<Matrix>::Ptr leftVector) :
Task(transformer.mName + ".MnaPostStep"), mTransformer(transformer), mLeftVector(leftVector) {
mAttributeDependencies.push_back(transformer.mSubInductor->attribute("i_intf"));
mAttributeDependencies.push_back(leftVector);
mModifiedAttributes.push_back(transformer.attribute("i_intf"));
mModifiedAttributes.push_back(transformer.attribute("v_intf"));
mTransformer.mnaAddPostStepDependencies(mPrevStepDependencies, mAttributeDependencies, mModifiedAttributes, mLeftVector);
}
void execute(Real time, Int timeStepCount);
void execute(Real time, Int timeStepCount) { mTransformer.mnaPostStep(time, timeStepCount, mLeftVector); };
private:
Transformer& mTransformer;
......
......@@ -11,8 +11,12 @@
using namespace CPS;
// #### General ####
SP::Ph1::Transformer::Transformer(String uid, String name, Logger::Level logLevel)
SP::Ph1::Transformer::Transformer(String uid, String name, Logger::Level logLevel, Bool withResistiveLosses)
: SimPowerComp<Complex>(uid, name, logLevel) {
if (withResistiveLosses)
setVirtualNodeNumber(3);
else
setVirtualNodeNumber(2);
mSLog->info("Create {} {}", this->type(), name);
mIntfVoltage = MatrixComp::Zero(1, 1);
......@@ -54,11 +58,6 @@ void SP::Ph1::Transformer::setParameters(Real nomVoltageEnd1, Real nomVoltageEnd
mNominalOmega = omega;
mReactance = mNominalOmega * mInductance;
if (mResistance > 0)
setVirtualNodeNumber(3);
else
setVirtualNodeNumber(2);
mSLog->info("Resistance={} [Ohm] Reactance={} [Ohm] (referred to primary side)", mResistance, mReactance);
mSLog->info("Tap Ratio={} [ ] Phase Shift={} [deg]", mRatioAbs, mRatioPhase);
......@@ -73,15 +72,81 @@ SimPowerComp<Complex>::Ptr SP::Ph1::Transformer::clone(String name) {
return copy;
}
void SP::Ph1::Transformer::initializeFromNodesAndTerminals(Real frequency) {
// Component parameters are referred to high voltage side.
// Switch terminals if transformer is connected the other way around.
if (Math::abs(mRatio) < 1.) {
mRatio = 1. / mRatio;
std::shared_ptr<SimTerminal<Complex>> tmp = mTerminals[0];
mTerminals[0] = mTerminals[1];
mTerminals[1] = tmp;
}
// Set initial voltage of virtual node in between
mVirtualNodes[0]->setInitialVoltage(initialSingleVoltage(1) * mRatio);
// Static calculations from load flow data
Real omega = 2. * PI * frequency;
Complex impedance = { mResistance, omega * mInductance };
mSLog->info("Reactance={} [Ohm] (referred to primary side)", omega * mInductance );
mIntfVoltage(0, 0) = mVirtualNodes[0]->initialSingleVoltage() - initialSingleVoltage(0);
mIntfCurrent(0, 0) = mIntfVoltage(0, 0) / impedance;
// Create series sub components
mSubInductor = std::make_shared<SP::Ph1::Inductor>(mUID + "_ind", mName + "_ind", Logger::Level::off);
mSubInductor->setParameters(mInductance);
if (mNumVirtualNodes == 3) {
mVirtualNodes[2]->setInitialVoltage(initialSingleVoltage(0));
mSubResistor = std::make_shared<SP::Ph1::Resistor>(mUID + "_res", mName + "_res", Logger::Level::off);
mSubResistor->setParameters(mResistance);
mSubResistor->connect({ node(0), mVirtualNodes[2] });
mSubResistor->initialize(mFrequencies);
mSubResistor->initializeFromNodesAndTerminals(frequency);
mSubInductor->connect({ mVirtualNodes[2], mVirtualNodes[0] });
} else {
mSubInductor->connect({ node(0), mVirtualNodes[0] });
}
mSubInductor->initialize(mFrequencies);
mSubInductor->initializeFromNodesAndTerminals(frequency);
// Create parallel sub components
// A snubber conductance is added on the low voltage side (resistance approximately scaled with LV side voltage)
if (Math::abs(mRatio) < 1.)
mSnubberResistance = std::abs(mNominalVoltageEnd1)*1e6;
else
mSnubberResistance = std::abs(mNominalVoltageEnd2)*1e6;
mSubSnubResistor = std::make_shared<SP::Ph1::Resistor>(mUID + "_snub_res", mName + "_snub_res", Logger::Level::off);
mSubSnubResistor->setParameters(mSnubberResistance);
mSubSnubResistor->connect({ node(1), SP::SimNode::GND });
mSubSnubResistor->initialize(mFrequencies);
mSubSnubResistor->initializeFromNodesAndTerminals(frequency);
mSLog->info("Snubber Resistance={} [Ohm] (connected to LV side)", mSnubberResistance);
mSLog->info(
"\n--- Initialization from powerflow ---"
"\nVoltage across: {:s}"
"\nCurrent: {:s}"
"\nTerminal 0 voltage: {:s}"
"\nTerminal 1 voltage: {:s}"
"\nVirtual Node 1 voltage: {:s}"
"\n--- Initialization from powerflow finished ---",
Logger::phasorToString(mIntfVoltage(0, 0)),
Logger::phasorToString(mIntfCurrent(0, 0)),
Logger::phasorToString(initialSingleVoltage(0)),
Logger::phasorToString(initialSingleVoltage(1)),
Logger::phasorToString(mVirtualNodes[0]->initialSingleVoltage()));
}
// #### Powerflow section ####
void SP::Ph1::Transformer::setBaseVoltage(Real baseVoltage) {
// Note: to be consistent set base voltage to higher voltage (and impedance values must be referred to high voltage side)
// TODO: use attribute setter for setting base voltage
mBaseVoltage = baseVoltage;
}
void SP::Ph1::Transformer::calculatePerUnitParameters(Real baseApparentPower, Real baseOmega) {
mSLog->info("#### Calculate Per Unit Parameters for {}", mName);
mBaseApparentPower = baseApparentPower;
......@@ -169,75 +234,11 @@ void SP::Ph1::Transformer::storeNodalInjection(Complex powerInjection) {
mStoreNodalPowerInjection = true;
}
// #### Getter ####
MatrixComp SP::Ph1::Transformer::Y_element() {
return mY_element;
}
// #### MNA Section ####
void SP::Ph1::Transformer::initializeFromNodesAndTerminals(Real frequency) {
// A snubber conductance is added on the low voltage side
Real snubberResistance = 1e3;
// Component parameters are referred to high voltage side.
// Switch terminals if transformer is connected the other way around.
if (Math::abs(mRatio) < 1.) {
mRatio = 1. / mRatio;
std::shared_ptr<SimTerminal<Complex>> tmp = mTerminals[0];
mTerminals[0] = mTerminals[1];
mTerminals[1] = tmp;
}
// Set initial voltage of virtual node in between
mVirtualNodes[0]->setInitialVoltage(initialSingleVoltage(1) * mRatio);
// Static calculations from load flow data
Real omega = 2. * PI * frequency;
Complex impedance = { mResistance, omega * mInductance };
mIntfVoltage(0, 0) = mVirtualNodes[0]->initialSingleVoltage() - initialSingleVoltage(0);
mIntfCurrent(0, 0) = mIntfVoltage(0, 0) / impedance;
// Create series sub components
mSubInductor = std::make_shared<SP::Ph1::Inductor>(mUID + "_ind", mName + "_ind", Logger::Level::off);
mSubInductor->setParameters(mInductance);
if (mNumVirtualNodes == 3) {
mVirtualNodes[2]->setInitialVoltage(initialSingleVoltage(0));
mSubResistor = std::make_shared<SP::Ph1::Resistor>(mUID + "_res", mName + "_res", Logger::Level::off);
mSubResistor->setParameters(mResistance);
mSubResistor->connect({ node(0), mVirtualNodes[2] });
mSubResistor->initializeFromNodesAndTerminals(frequency);
mSubInductor->connect({ mVirtualNodes[2], mVirtualNodes[0] });
}
else {
mSubInductor->connect({ node(0), mVirtualNodes[0] });
}
mSubInductor->initializeFromNodesAndTerminals(frequency);
// Create parallel sub components
mSubSnubResistor = std::make_shared<SP::Ph1::Resistor>(mUID + "_snub_res", mName + "_snub_res", Logger::Level::off);
mSubSnubResistor->setParameters(snubberResistance);
mSubSnubResistor->connect({ node(1), SP::SimNode::GND });
mSubSnubResistor->initializeFromNodesAndTerminals(frequency);
mSLog->info(
"\n--- Initialization from powerflow ---"
"\nVoltage across: {:s}"
"\nCurrent: {:s}"
"\nTerminal 0 voltage: {:s}"
"\nTerminal 1 voltage: {:s}"
"\nVirtual Node 1 voltage: {:s}"
"\n--- Initialization from powerflow finished ---",
Logger::phasorToString(mIntfVoltage(0, 0)),
Logger::phasorToString(mIntfCurrent(0, 0)),
Logger::phasorToString(initialSingleVoltage(0)),
Logger::phasorToString(initialSingleVoltage(1)),
Logger::phasorToString(mVirtualNodes[0]->initialSingleVoltage()));
}
void SP::Ph1::Transformer::mnaInitialize(Real omega, Real timeStep, Attribute<Matrix>::Ptr leftVector) {
MNAInterface::mnaInitialize(omega, timeStep);
......@@ -245,13 +246,16 @@ void SP::Ph1::Transformer::mnaInitialize(Real omega, Real timeStep, Attribute<Ma
mRightVector = Matrix::Zero(leftVector->get().rows(), 1);
auto subComponents = MNAInterface::List({ mSubInductor, mSubSnubResistor });
if (mSubResistor)
subComponents.push_back(mSubResistor);
for (auto comp : subComponents) {
comp->mnaInitialize(omega, timeStep, leftVector);
for (auto task : comp->mnaTasks()) {
for (auto task : comp->mnaTasks())
mMnaTasks.push_back(task);
}
}
mMnaTasks.push_back(std::make_shared<MnaPostStep>(*this, leftVector));
......@@ -297,11 +301,28 @@ void SP::Ph1::Transformer::mnaApplySystemMatrixStamp(Matrix& systemMatrix) {
}
void SP::Ph1::Transformer::MnaPostStep::execute(Real time, Int timeStepCount) {
mTransformer.mnaUpdateVoltage(*mLeftVector);
mTransformer.mnaUpdateCurrent(*mLeftVector);
void SP::Ph1::Transformer::mnaAddPostStepDependencies(AttributeBase::List &prevStepDependencies, AttributeBase::List &attributeDependencies, AttributeBase::List &modifiedAttributes, Attribute<Matrix>::Ptr &leftVector) {
// add post-step dependencies of subcomponents
if (mSubResistor)
this->mSubResistor->mnaAddPostStepDependencies(prevStepDependencies, attributeDependencies, modifiedAttributes, leftVector);
this->mSubInductor->mnaAddPostStepDependencies(prevStepDependencies, attributeDependencies, modifiedAttributes, leftVector);
this->mSubSnubResistor->mnaAddPostStepDependencies(prevStepDependencies, attributeDependencies, modifiedAttributes, leftVector);
// add post-step dependencies of component itself
attributeDependencies.push_back(leftVector);
modifiedAttributes.push_back(this->attribute("v_intf"));
modifiedAttributes.push_back(this->attribute("i_intf"));
}
void SP::Ph1::Transformer::mnaPostStep(Real time, Int timeStepCount, Attribute<Matrix>::Ptr &leftVector) {
// post-step of subcomponents
if (mSubResistor)
this->mSubResistor->mnaPostStep(time, timeStepCount, leftVector);
this->mSubInductor->mnaPostStep(time, timeStepCount, leftVector);
this->mSubSnubResistor->mnaPostStep(time, timeStepCount, leftVector);
// post-step of component itself
this->mnaUpdateVoltage(*leftVector);
this->mnaUpdateCurrent(*leftVector);
}
void SP::Ph1::Transformer::mnaUpdateCurrent(const Matrix& leftVector) {
mIntfCurrent(0, 0) = mSubInductor->intfCurrent()(0, 0);
......@@ -309,14 +330,10 @@ void SP::Ph1::Transformer::mnaUpdateCurrent(const Matrix& leftVector) {
}
void SP::Ph1::Transformer::mnaUpdateVoltage(const Matrix& leftVector) {
// TODO is this correct?
mIntfVoltage(0, 0) = mVirtualNodes[0]->initialSingleVoltage() - initialSingleVoltage(0);
// v1 - v0
mIntfVoltage(0, 0) = 0;
mIntfVoltage(0, 0) = Math::complexFromVectorElement(leftVector, matrixNodeIndex(1));
mIntfVoltage(0, 0) = mIntfVoltage(0, 0) - Math::complexFromVectorElement(leftVector, mVirtualNodes[0]->matrixNodeIndex());
mSLog->debug("Voltage {:s}", Logger::phasorToString(mIntfVoltage(0, 0)));
}
Markdown is supported
0%