Merge branch 'decoupling-line' into 'master'

Decoupling line

See merge request acs/public/simulation/dpsim!130

Former-commit-id: 97ce49fa

Dependencies/libcps/Include/cps/Attribute.h

-/** Attributes
- *
- * @file
- * @author Steffen Vogel <stvogel@eonerc.rwth-aachen.de>
- * @copyright 2018, Institute for Automation of Complex Power Systems, EONERC
+/**
+ * @copyright 2017 Institute for Automation of Complex Power Systems, EONERC
  *
  * CPowerSystems
  *

Dependencies/libcps/Include/cps/Components.h

@@ -112,6 +112,7 @@
 #include <cps/EMT/EMT_Ph3_Switch.h>
 
 #include <cps/Signal/DecouplingLine.h>
+#include <cps/Signal/DecouplingLineEMT.h>
 #include <cps/Signal/Exciter.h>
 #include <cps/Signal/TurbineGovernor.h>
 #include <cps/Signal/FIRFilter.h>

Dependencies/libcps/Include/cps/Signal/DecouplingLine.h

@@ -43,7 +43,7 @@ namespace Signal {
 		Complex mSrcCur1Ref;
 		Complex mSrcCur2Ref;
 
-		// TODO also EMT support?
+		std::shared_ptr<DP::Node> mNode1, mNode2;
 		std::shared_ptr<DP::Ph1::Resistor> mRes1, mRes2;
 		std::shared_ptr<DP::Ph1::CurrentSource> mSrc1, mSrc2;
 		Attribute<Complex>::Ptr mSrcCur1, mSrcCur2;
@@ -64,6 +64,9 @@ namespace Signal {
 		DecouplingLine(String name, Node<Complex>::Ptr node1, Node<Complex>::Ptr node2,
 			Real resistance, Real inductance, Real capacitance, Logger::Level logLevel = Logger::Level::info);
 
+		DecouplingLine(String name, Logger::Level logLevel = Logger::Level::info);
+
+		void setParameters(Node<Complex>::Ptr node1, Node<Complex>::Ptr node2, Real resistance, Real inductance, Real capacitance);
 		void initialize(Real omega, Real timeStep);
 		void step(Real time, Int timeStepCount);
 		void postStep();

Dependencies/libcps/Include/cps/Signal/DecouplingLineEMT.h

+/**
+ * @copyright 2017 Institute for Automation of Complex Power Systems, EONERC
+ *
+ * CPowerSystems
+ *
+ * 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/>.
+ *********************************************************************************/
+
+#pragma once
+
+#include <vector>
+
+#include <cps/EMT/EMT_Ph1_CurrentSource.h>
+#include <cps/EMT/EMT_Ph1_Resistor.h>
+#include <cps/SignalComponent.h>
+#include <cps/Task.h>
+
+namespace CPS {
+namespace Signal {
+	class DecouplingLineEMT :
+		public SignalComponent,
+		public SharedFactory<DecouplingLineEMT> {
+	protected:
+		Real mDelay;
+		Real mResistance;
+		Real mInductance;
+		Real mCapacitance;
+		Real mSurgeImpedance;
+		Real mSrcCur1Ref;
+		Real mSrcCur2Ref;
+
+		std::shared_ptr<EMT::Node> mNode1, mNode2;
+		std::shared_ptr<EMT::Ph1::Resistor> mRes1, mRes2;
+		std::shared_ptr<EMT::Ph1::CurrentSource> mSrc1, mSrc2;
+		Attribute<Complex>::Ptr mSrcCur1, mSrcCur2;
+
+		// Ringbuffers for the values of previous timesteps
+		// TODO make these matrix attributes
+		std::vector<Real> mVolt1, mVolt2, mCur1, mCur2;
+		// workaround for dependency analysis as long as the states aren't attributes
+		Matrix mStates;
+		UInt mBufIdx = 0;
+		UInt mBufSize;
+		Real mAlpha;
+
+		Real interpolate(std::vector<Real>& data);
+	public:
+		typedef std::shared_ptr<DecouplingLineEMT> Ptr;
+
+		DecouplingLineEMT(String name, Logger::Level logLevel = Logger::Level::info);
+
+		void setParameters(Node<Real>::Ptr node1, Node<Real>::Ptr node2,
+			Real resistance, Real inductance, Real capacitance);
+		void initialize(Real omega, Real timeStep);
+		void step(Real time, Int timeStepCount);
+		void postStep();
+		Task::List getTasks();
+		Component::List getLineComponents();
+
+		class PreStep : public Task {
+		public:
+			PreStep(DecouplingLineEMT& line) :
+				Task(line.mName + ".MnaPreStep"), mLine(line) {
+				mPrevStepDependencies.push_back(mLine.attribute("states"));
+				mModifiedAttributes.push_back(mLine.mSrc1->attribute("I_ref"));
+				mModifiedAttributes.push_back(mLine.mSrc2->attribute("I_ref"));
+			}
+
+			void execute(Real time, Int timeStepCount);
+
+		private:
+			DecouplingLineEMT& mLine;
+		};
+
+		class PostStep : public Task {
+		public:
+			PostStep(DecouplingLineEMT& line) :
+				Task(line.mName + ".PostStep"), mLine(line) {
+				mAttributeDependencies.push_back(mLine.mRes1->attribute("v_intf"));
+				mAttributeDependencies.push_back(mLine.mRes1->attribute("i_intf"));
+				mAttributeDependencies.push_back(mLine.mRes2->attribute("v_intf"));
+				mAttributeDependencies.push_back(mLine.mRes2->attribute("i_intf"));
+				mModifiedAttributes.push_back(mLine.attribute("states"));
+			}
+
+			void execute(Real time, Int timeStepCount);
+
+		private:
+			DecouplingLineEMT& mLine;
+		};
+	};
+}
+}

Dependencies/libcps/Source/CMakeLists.txt

@@ -111,8 +111,8 @@ list(APPEND SOURCES
 	SP/SP_Ph3_VoltageSource.cpp
 	SP/SP_Ph3_ControlledVoltageSource.cpp
 
-
 	Signal/DecouplingLine.cpp
+	Signal/DecouplingLineEMT.cpp
 	Signal/Exciter.cpp
 	Signal/FIRFilter.cpp
 	Signal/TurbineGovernor.cpp

Dependencies/libcps/Source/Logger.cpp

 /**
- * @author Markus Mirz <mmirz@eonerc.rwth-aachen.de>
- * @copyright 2017-2018, Institute for Automation of Complex Power Systems, EONERC
+ * @copyright 2017 Institute for Automation of Complex Power Systems, EONERC
  *
  * CPowerSystems
  *

Dependencies/libcps/Source/Signal/DecouplingLine.cpp

-/** Signalling part of a decoupling transmission line
- *
- * @file
- * @author Georg Reinke <georg.reinke@rwth-aachen.de>
- * @copyright 2017-2018, Institute for Automation of Complex Power Systems, EONERC
+/**
+ * @copyright 2017 Institute for Automation of Complex Power Systems, EONERC
  *
  * CPowerSystems
  *
@@ -30,6 +27,7 @@ DecouplingLine::DecouplingLine(String name, Node<Complex>::Ptr node1, Node<Compl
 	Real resistance, Real inductance, Real capacitance, Logger::Level logLevel) :
 	SignalComponent(name, name, logLevel),
 	mResistance(resistance), mInductance(inductance), mCapacitance(capacitance) {
+
 	addAttribute<Matrix>("states", &mStates);
 	addAttribute<Complex>("i_src1", &mSrcCur1Ref, Flags::read);
 	addAttribute<Complex>("i_src2", &mSrcCur2Ref, Flags::read);
@@ -56,6 +54,46 @@ DecouplingLine::DecouplingLine(String name, Node<Complex>::Ptr node1, Node<Compl
 	mSrcCur2 = mSrc2->attributeComplex("I_ref");
 }
 
+DecouplingLine::DecouplingLine(String name, Logger::Level logLevel) :
+	SignalComponent(name, name, logLevel) {
+
+	addAttribute<Matrix>("states", &mStates);
+	addAttribute<Complex>("i_src1", &mSrcCur1Ref, Flags::read);
+	addAttribute<Complex>("i_src2", &mSrcCur2Ref, Flags::read);
+
+	mRes1 = Resistor::make(name + "_r1", logLevel);
+	mRes2 = Resistor::make(name + "_r2", logLevel);
+	mSrc1 = CurrentSource::make(name + "_i1", logLevel);
+	mSrc2 = CurrentSource::make(name + "_i2", logLevel);
+
+	mSrcCur1 = mSrc1->attributeComplex("I_ref");
+	mSrcCur2 = mSrc2->attributeComplex("I_ref");
+}
+
+void DecouplingLine::setParameters(Node<Complex>::Ptr node1, Node<Complex>::Ptr node2,
+	Real resistance, Real inductance, Real capacitance) {
+
+	mResistance = resistance;
+	mInductance = inductance;
+	mCapacitance = capacitance;
+	mNode1 = node1;
+	mNode2 = node2;
+
+	mSurgeImpedance = sqrt(inductance / capacitance);
+	mDelay = sqrt(inductance * capacitance);
+	mSLog->info("surge impedance: {}", mSurgeImpedance);
+	mSLog->info("delay: {}", mDelay);
+
+	mRes1->setParameters(mSurgeImpedance + resistance / 4);
+	mRes1->connect({node1, Node<Complex>::GND});
+	mRes2->setParameters(mSurgeImpedance + resistance / 4);
+	mRes2->connect({node2, Node<Complex>::GND});
+	mSrc1->setParameters(0);
+	mSrc1->connect({node1, Node<Complex>::GND});
+	mSrc2->setParameters(0);
+	mSrc2->connect({node2, Node<Complex>::GND});
+}
+
 void DecouplingLine::initialize(Real omega, Real timeStep) {
 	if (mDelay < timeStep)
 		throw SystemError("Timestep too large for decoupling");
@@ -64,14 +102,16 @@ void DecouplingLine::initialize(Real omega, Real timeStep) {
 	mAlpha = 1 - (mBufSize - mDelay / timeStep);
 	mSLog->info("bufsize {} alpha {}", mBufSize, mAlpha);
 
-	Complex volt1 = mRes1->initialSingleVoltage(0);
-	Complex volt2 = mRes2->initialSingleVoltage(0);
+	Complex volt1 = mNode1->initialSingleVoltage();
+	Complex volt2 = mNode2->initialSingleVoltage();
 	// TODO different initialization for lumped resistance?
 	Complex initAdmittance = 1. / Complex(mResistance, omega * mInductance) + Complex(0, omega * mCapacitance / 2);
 	Complex cur1 = volt1 * initAdmittance - volt2 / Complex(mResistance, omega * mInductance);
 	Complex cur2 = volt2 * initAdmittance - volt1 / Complex(mResistance, omega * mInductance);
 	mSLog->info("initial voltages: v_k {} v_m {}", volt1, volt2);
 	mSLog->info("initial currents: i_km {} i_mk {}", cur1, cur2);
+
+	// Resize ring buffers and initialize
 	mVolt1.resize(mBufSize, volt1);
 	mVolt2.resize(mBufSize, volt2);
 	mCur1.resize(mBufSize, cur1);
@@ -94,17 +134,19 @@ void DecouplingLine::step(Real time, Int timeStepCount) {
 	if (timeStepCount == 0) {
 		// bit of a hack for proper initialization
 		mSrcCur1Ref = cur1 - volt1 / (mSurgeImpedance + mResistance / 4);
-		mSrcCur1Ref = cur2 - volt2 / (mSurgeImpedance + mResistance / 4);
+		mSrcCur2Ref = cur2 - volt2 / (mSurgeImpedance + mResistance / 4);
 	} else {
 		// Update currents
 		Real denom = (mSurgeImpedance + mResistance/4) * (mSurgeImpedance + mResistance/4);
 		mSrcCur1Ref = -mSurgeImpedance / denom * (volt2 + (mSurgeImpedance - mResistance/4) * cur2)
 			- mResistance/4 / denom * (volt1 + (mSurgeImpedance - mResistance/4) * cur1);
-		mSrcCur1Ref = -mSurgeImpedance / denom * (volt1 + (mSurgeImpedance - mResistance/4) * cur1)
+		mSrcCur2Ref = -mSurgeImpedance / denom * (volt1 + (mSurgeImpedance - mResistance/4) * cur1)
 			- mResistance/4 / denom * (volt2 + (mSurgeImpedance - mResistance/4) * cur2);
+		mSrcCur1Ref = mSrcCur1Ref * Complex(cos(-2.*PI*50*mDelay),sin(-2.*PI*50*mDelay));
+		mSrcCur2Ref = mSrcCur2Ref * Complex(cos(-2.*PI*50*mDelay),sin(-2.*PI*50*mDelay));
 	}
 	mSrcCur1->set(mSrcCur1Ref);
-	mSrcCur2->set(mSrcCur1Ref);
+	mSrcCur2->set(mSrcCur2Ref);
 }
 
 void DecouplingLine::PreStep::execute(Real time, Int timeStepCount) {

Dependencies/libcps/Source/Signal/DecouplingLineEMT.cpp

+/**
+ * @copyright 2017 Institute for Automation of Complex Power Systems, EONERC
+ *
+ * CPowerSystems
+ *
+ * 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 <cps/Signal/DecouplingLineEMT.h>
+
+using namespace CPS;
+using namespace CPS::EMT::Ph1;
+using namespace CPS::Signal;
+
+DecouplingLineEMT::DecouplingLineEMT(String name, Logger::Level logLevel) :
+	SignalComponent(name, name, logLevel) {
+
+	addAttribute<Matrix>("states", &mStates);
+	addAttribute<Real>("i_src1", &mSrcCur1Ref, Flags::read);
+	addAttribute<Real>("i_src2", &mSrcCur2Ref, Flags::read);
+
+	mRes1 = Resistor::make(name + "_r1", logLevel);
+	mRes2 = Resistor::make(name + "_r2", logLevel);
+	mSrc1 = CurrentSource::make(name + "_i1", logLevel);
+	mSrc2 = CurrentSource::make(name + "_i2", logLevel);
+
+	mSrcCur1 = mSrc1->attributeComplex("I_ref");
+	mSrcCur2 = mSrc2->attributeComplex("I_ref");
+}
+
+void DecouplingLineEMT::setParameters(Node<Real>::Ptr node1, Node<Real>::Ptr node2,
+	Real resistance, Real inductance, Real capacitance) {
+
+	mResistance = resistance;
+	mInductance = inductance;
+	mCapacitance = capacitance;
+	mNode1 = node1;
+	mNode2 = node2;
+
+	mSurgeImpedance = sqrt(inductance / capacitance);
+	mDelay = sqrt(inductance * capacitance);
+	mSLog->info("surge impedance: {}", mSurgeImpedance);
+	mSLog->info("delay: {}", mDelay);
+
+	mRes1->setParameters(mSurgeImpedance + mResistance / 4);
+	mRes1->connect({node1, Node<Real>::GND});
+	mRes2->setParameters(mSurgeImpedance + mResistance / 4);
+	mRes2->connect({node2, Node<Real>::GND});
+	mSrc1->setParameters(0);
+	mSrc1->connect({node1, Node<Real>::GND});
+	mSrc2->setParameters(0);
+	mSrc2->connect({node2, Node<Real>::GND});
+}
+
+void DecouplingLineEMT::initialize(Real omega, Real timeStep) {
+	if (mDelay < timeStep)
+		throw SystemError("Timestep too large for decoupling");
+
+	mBufSize = static_cast<UInt>(ceil(mDelay / timeStep));
+	mAlpha = 1 - (mBufSize - mDelay / timeStep);
+	mSLog->info("bufsize {} alpha {}", mBufSize, mAlpha);
+
+	// Initialization based on static PI-line model
+	Complex volt1 = mNode1->initialSingleVoltage();
+	Complex volt2 = mNode2->initialSingleVoltage();
+	Complex initAdmittance = 1. / Complex(mResistance, omega * mInductance) + Complex(0, omega * mCapacitance / 2);
+	Complex cur1 = volt1 * initAdmittance - volt2 / Complex(mResistance, omega * mInductance);
+	Complex cur2 = volt2 * initAdmittance - volt1 / Complex(mResistance, omega * mInductance);
+	mSLog->info("initial voltages: v_k {} v_m {}", volt1, volt2);
+	mSLog->info("initial currents: i_km {} i_mk {}", cur1, cur2);
+
+	// Resize ring buffers and initialize
+	mVolt1.resize(mBufSize, volt1.real());
+	mVolt2.resize(mBufSize, volt2.real());
+	mCur1.resize(mBufSize, cur1.real());
+	mCur2.resize(mBufSize, cur2.real());
+}
+
+Real DecouplingLineEMT::interpolate(std::vector<Real>& data) {
+	// linear interpolation of the nearest values
+	Real c1 = data[mBufIdx];
+	Real c2 = mBufIdx == mBufSize-1 ? data[0] : data[mBufIdx+1];
+	return mAlpha * c1 + (1-mAlpha) * c2;
+}
+
+void DecouplingLineEMT::step(Real time, Int timeStepCount) {
+	Real volt1 = interpolate(mVolt1);
+	Real volt2 = interpolate(mVolt2);
+	Real cur1 = interpolate(mCur1);
+	Real cur2 = interpolate(mCur2);
+	Real denom = (mSurgeImpedance + mResistance/4) * (mSurgeImpedance + mResistance/4);
+
+	if (timeStepCount == 0) {
+		// initialization
+		mSrcCur1Ref = cur1 - volt1 / (mSurgeImpedance + mResistance / 4);
+		mSrcCur2Ref = cur2 - volt2 / (mSurgeImpedance + mResistance / 4);
+	} else {
+		// Update currents
+		mSrcCur1Ref = -mSurgeImpedance / denom * (volt2 + (mSurgeImpedance - mResistance/4) * cur2)
+			-mResistance/4 / denom * (volt1 + (mSurgeImpedance - mResistance/4) * cur1);
+		mSrcCur2Ref = -mSurgeImpedance / denom * (volt1 + (mSurgeImpedance - mResistance/4) * cur1)
+			-mResistance/4 / denom * (volt2 + (mSurgeImpedance - mResistance/4) * cur2);
+	}
+	mSrcCur1->set(mSrcCur1Ref);
+	mSrcCur2->set(mSrcCur2Ref);
+}
+
+void DecouplingLineEMT::PreStep::execute(Real time, Int timeStepCount) {
+	mLine.step(time, timeStepCount);
+}
+
+void DecouplingLineEMT::postStep() {
+	// Update ringbuffers with new values
+	mVolt1[mBufIdx] = -mRes1->intfVoltage()(0,0);
+	mVolt2[mBufIdx] = -mRes2->intfVoltage()(0,0);
+	mCur1[mBufIdx] = -mRes1->intfCurrent()(0,0) + mSrcCur1->get().real();
+	mCur2[mBufIdx] = -mRes2->intfCurrent()(0,0) + mSrcCur2->get().real();
+
+	mBufIdx++;
+	if (mBufIdx == mBufSize)
+		mBufIdx = 0;
+}
+
+void DecouplingLineEMT::PostStep::execute(Real time, Int timeStepCount) {
+	mLine.postStep();
+}
+
+Task::List DecouplingLineEMT::getTasks() {
+	return Task::List({std::make_shared<PreStep>(*this), std::make_shared<PostStep>(*this)});
+}
+
+Component::List DecouplingLineEMT::getLineComponents() {
+	return Component::List({mRes1, mRes2, mSrc1, mSrc2});
+}

Examples/Cxx/Circuits/DP_DecouplingLine.cpp

 /**
- * @file
- * @author Markus Mirz <mmirz@eonerc.rwth-aachen.de>
- * @copyright 2017-2019, Institute for Automation of Complex Power Systems, EONERC
+ * @copyright 2017 Institute for Automation of Complex Power Systems, EONERC
  *
  * DPsim
  *
@@ -92,15 +90,14 @@ void simDecoupling() {
 	auto n2 = Node::make("n2");
 
 	// Components
-	auto vs = Ph1::VoltageSource::make("v_1");
-	vs->setParameters(CPS::Math::polar(100000, -PI/2.));
+	auto vs = Ph1::VoltageSource::make("Vsrc");
+	vs->setParameters(CPS::Math::polar(100000, 0));
 
 	Real resistance = 5;
 	Real inductance = 0.16;
 	Real capacitance = 1.0e-6;
-	auto dline = CPS::Signal::DecouplingLine::make("dline",
-		n1, n2, resistance, inductance, capacitance,
-		Logger::Level::debug);
+	auto dline = CPS::Signal::DecouplingLine::make("DecLine", Logger::Level::debug);
+	dline->setParameters(n1, n2, resistance, inductance, capacitance);
 
 	auto load = Ph1::Resistor::make("R_load");
 	load->setParameters(10000);
@@ -133,7 +130,6 @@ void simDecoupling() {
 }
 
 void simDecouplingEMT() {
-	/*
 	Real timeStep = 0.00005;
 	Real finalTime = 0.1;
 	String simName = "EMT_Decoupling_Wave";
@@ -144,16 +140,14 @@ void simDecouplingEMT() {
 	auto n2 = CPS::EMT::Node::make("n2");
 
 	// Components
-	auto vs = CPS::EMT::Ph1::VoltageSource::make("v_1");
-	vs->setParameters(CPS::Math::polar(100000, -PI/2.));
-
+	auto vs = CPS::EMT::Ph1::VoltageSource::make("Vsrc_emt");
+	vs->setParameters(CPS::Math::polar(100000, 0), 50);
 
 	Real resistance = 5;
 	Real inductance = 0.16;
 	Real capacitance = 1.0e-6;
-	auto dline = CPS::Signal::DecouplingLine::make("dline",
-		n1, n2, resistance, inductance, capacitance,
-		Logger::Level::debug);
+	auto dline = CPS::Signal::DecouplingLineEMT::make("DecLine_emt", Logger::Level::debug);
+	dline->setParameters(n1, n2, resistance, inductance, capacitance);
 
 	auto load = CPS::EMT::Ph1::Resistor::make("R_load");
 	load->setParameters(10000);
@@ -173,6 +167,8 @@ void simDecouplingEMT() {
 	logger->addAttribute("v2", n2->attribute("v"));
 	logger->addAttribute("i1", vs->attribute("i_intf"));
 	logger->addAttribute("i2", load->attribute("i_intf"));
+	logger->addAttribute("i_src1", dline->attribute("i_src1"));
+	logger->addAttribute("i_src2", dline->attribute("i_src2"));
 
 	Simulation sim(simName);
 	sim.setSystem(sys);
@@ -182,7 +178,6 @@ void simDecouplingEMT() {
 	sim.addLogger(logger);
 
 	sim.run();
-	*/
 }
 
 int main(int argc, char* argv[]) {

Examples/Cxx/Circuits/DP_PiLine.cpp

-/** Reference Circuits
- *
- * @author Markus Mirz <mmirz@eonerc.rwth-aachen.de>
- * @copyright 2017-2018, Institute for Automation of Complex Power Systems, EONERC
+/**
+ * @copyright 2017 Institute for Automation of Complex Power Systems, EONERC
  *
  * DPsim
  *

Examples/Notebooks/Components/Line.ipynb.REMOVED.git-id

-7fee9a85f6541ec5b1586d7265b0f366e24f1a82
\ No newline at end of file
+e7495b8f7843e57c3e3cd8ece9bee7e00b1cfd1a
\ No newline at end of file