PMU signal implementation using GSL

dpsim-models/include/dpsim-models/Components.h

@@ -131,3 +131,4 @@
 #include <dpsim-models/Signal/SineWaveGenerator.h>
 #include <dpsim-models/Signal/FrequencyRampGenerator.h>
 #include <dpsim-models/Signal/CosineFMGenerator.h>
+#include <dpsim-models/Signal/PMUSignalDevice.h>

dpsim-models/include/dpsim-models/Signal/PMUSignalDevice.h

@@ -0,0 +1,55 @@
+
+#pragma once
+
+#include <vector>
+
+#include <dpsim-models/SimPowerComp.h>
+#include <dpsim-models/SimSignalComp.h>
+#include <dpsim-models/Task.h>
+
+
+namespace CPS {
+namespace Signal {
+
+    class PMUSignalDevice:
+        public SimSignalComp,
+		public SharedFactory<PMUSignalDevice> 
+        {
+	protected:
+	public:
+		///FIXME: This is never explicitely set to reference anything, so the outside code is responsible for setting up the reference.
+
+		//***PMUSignalDevice has two attributes:Input and Output***//
+		const  Attribute<MatrixComp>::Ptr mInput;
+        const  Attribute<MatrixComp>::Ptr mOutput;
+		Real mSigma;
+
+		//Constructor 
+		PMUSignalDevice(String name, Logger::Level logLevel = Logger::Level::off);
+
+		//Operation for adding measurement error
+		void MeasurementError(Real time);
+
+		//Set the gaussian standard deviation for this PMU
+	    void setParameters(Real Sigma);
+
+		//***Get task list for the solver***//
+		Task::List getTasks();
+
+		//***The process of the solve is recognized as a PostStep.***//
+ 		class PostStep : public Task {
+			public:
+				PostStep(PMUSignalDevice& PMU) :
+					Task(**PMU.mName + ".Poststep"), mPMU(PMU) {
+						mAttributeDependencies.push_back(PMU.mInput);
+						mModifiedAttributes.push_back(PMU.mOutput);
+					}
+
+				void execute(Real time, Int timeStepCount);
+
+			private:
+				PMUSignalDevice& mPMU;
+			};
+        };
+    }
+}

dpsim-models/src/CMakeLists.txt

@@ -146,6 +146,7 @@ list(APPEND MODELS_SOURCES
 	Signal/SignalGenerator.cpp
 	Signal/FrequencyRampGenerator.cpp
 	Signal/CosineFMGenerator.cpp
+	Signal/PMUSignalDevice.cpp
 )
 
 if(WITH_CIM)

dpsim-models/src/Signal/PMUSignalDevice.cpp

@@ -0,0 +1,57 @@
+#include <dpsim-models/Signal/PMUSignalDevice.h>
+#include <gsl/gsl_rng.h>
+#include <gsl/gsl_randist.h>
+#include <time.h>
+using namespace CPS;
+using namespace CPS::Signal;
+
+PMUSignalDevice::PMUSignalDevice(String name, Logger::Level logLevel) :
+	SimSignalComp(name, name, logLevel),
+    /// CHECK: Which of these really need to be attributes?
+    mInput(mAttributes->createDynamic<MatrixComp>("input")), 
+
+    //*** Important: the output attribute is a matrix component, which should be initialized as Zero(1,1)
+    mOutput(mAttributes->create<MatrixComp>("output", MatrixComp::Zero(1,1))){ }
+
+
+// }
+
+void PMUSignalDevice::setParameters(Real Sigma){
+    mSigma=Sigma;
+}
+
+void PMUSignalDevice::MeasurementError(Real time){
+
+    // Relative standard measurement uncertainty
+    double relMeasUncertainty;
+
+    // GSL random number generator
+    gsl_rng *r=gsl_rng_alloc(gsl_rng_default);
+
+    // To avoid the seed of random number always selected as the same, we use the UTC time to generate the seed.
+    // Time consumption of the calculations betweens elements is small, here the "tv_nsec" can distinct the time difference within nano second.
+    struct timespec ts;
+    timespec_get(&ts,TIME_UTC);
+    gsl_rng_set(r,ts.tv_nsec);
+
+    // The relative measurement uncertainty is gaussian distributed, mean value 0, standard deviation mSigma.
+    relMeasUncertainty = gsl_ran_gaussian (r, mSigma);
+
+    // Output = Input + error 
+    //        = Input + k*Sigma 
+    //        = Input + Input*relMeasUncertainty
+     **mOutput =**mInput*(1+relMeasUncertainty);
+
+    // To log the GSL random number, mInput and mOutput.
+    SPDLOG_LOGGER_INFO(mSLog, "GSL value: gsl= {}", relMeasUncertainty);
+    SPDLOG_LOGGER_INFO(mSLog,"Input values: input = {}", **mInput);
+    SPDLOG_LOGGER_INFO(mSLog,"Output values: output = {}:", **mOutput);
+}
+
+void PMUSignalDevice::PostStep::execute(Real time, Int timeStepCount) {
+ 	mPMU.MeasurementError(time);
+}
+
+Task::List PMUSignalDevice::getTasks() {
+	return Task::List({std::make_shared<PMUSignalDevice::PostStep>(*this)});
+}

dpsim/examples/cxx/CMakeLists.txt

@@ -89,6 +89,8 @@ set(CIRCUIT_SOURCES
 	Circuits/DP_SynGenTrStab_3Bus_Fault.cpp
 	Circuits/SP_SynGenTrStab_3Bus_Fault.cpp
 
+	Circuits/Three_bus_sim.cpp
+
 )
 
 set(SYNCGEN_SOURCES

dpsim/examples/cxx/Circuits/Three_bus_sim.cpp

@@ -0,0 +1,149 @@
+#include <DPsim.h>
+#include "../Examples.h"
+
+using namespace DPsim;
+using namespace CPS;
+using namespace CIM::Examples::Grids::Three_bus_cfg;
+
+ScenarioConfig Three_bus_cfg;
+
+void Three_bus_sim(String simName, Real timeStep, Real finalTime, Real cmdInertia_G1, Real cmdInertia_G2, Real cmdDamping_G1, Real cmdDamping_G2)
+{
+	// ----- POWERFLOW FOR INITIALIZATION -----
+	Real timeStepPF = finalTime;
+	Real finalTimePF = finalTime + timeStepPF;
+	String simNamePF = simName + "_PF";
+	Logger::setLogDir("logs/" + simNamePF);
+
+	// Components
+	// nodes
+	auto n1_PF = SimNode<Complex>::make("n1", PhaseType::Single);
+	auto n2_PF = SimNode<Complex>::make("n2", PhaseType::Single);
+	auto n3_PF = SimNode<Complex>::make("n3", PhaseType::Single);
+
+	//***Applying the new signal model: PMUSignalDevice
+	//***Tonstruction a PMU
+	auto pmu_1 = Signal::PMUSignalDevice::make("PMU_1", Logger::Level::debug);
+	//***The attribute mInput of PMU is set by the mVoltage of the SimNode
+	pmu_1->mInput->setReference(n1_PF->mVoltage);
+
+	//*** PMU signal deveice produces a Gaussian distributed measurement error rate***//
+	//*** K * Sigma = real value * errorRate                                       ***//
+	pmu_1->setParameters(0.01);
+
+	auto pmu_2 = Signal::PMUSignalDevice::make("PMU_2", Logger::Level::debug);
+	pmu_2->mInput->setReference(n2_PF->mVoltage);
+	pmu_2->setParameters(0.001);
+
+	auto pmu_3 = Signal::PMUSignalDevice::make("PMU_3", Logger::Level::debug);
+	pmu_3->mInput->setReference(n3_PF->mVoltage);
+	pmu_3->setParameters(0.001);
+
+	// injection
+	// configuration(Synchronous generator 1)
+	auto gen1_PF = SP::Ph1::SynchronGenerator::make("Generator_1", Logger::Level::debug);
+	// setPointVoltage is defined as the voltage at the transfomer primary side and should be transformed to network side
+	gen1_PF->setParameters(Three_bus_cfg.nomPower_G1, Three_bus_cfg.nomPhPhVoltRMS_G1, Three_bus_cfg.initActivePower_G1, Three_bus_cfg.setPointVoltage_G1 * Three_bus_cfg.t1_ratio, PowerflowBusType::VD, Three_bus_cfg.initReactivePower_G1);
+	gen1_PF->setBaseVoltage(Three_bus_cfg.Vnom);
+
+	// Synchronous generator 2
+	auto gen2_PF = SP::Ph1::SynchronGenerator::make("Generator_2", Logger::Level::debug);
+	// setPointVoltage is defined as the voltage at the transfomer primary side and should be transformed to network side
+	gen2_PF->setParameters(Three_bus_cfg.nomPower_G2, Three_bus_cfg.nomPhPhVoltRMS_G2, Three_bus_cfg.initActivePower_G2, Three_bus_cfg.setPointVoltage_G2 * Three_bus_cfg.t2_ratio, PowerflowBusType::PV, Three_bus_cfg.initReactivePower_G2);
+	gen2_PF->setBaseVoltage(Three_bus_cfg.Vnom);
+
+	// use Shunt as Load for powerflow
+	auto load_PF = SP::Ph1::Shunt::make("Load", Logger::Level::debug);
+	load_PF->setParameters(Three_bus_cfg.activePower_L / std::pow(Three_bus_cfg.Vnom, 2), -Three_bus_cfg.reactivePower_L / std::pow(Three_bus_cfg.Vnom, 2));
+	load_PF->setBaseVoltage(Three_bus_cfg.Vnom);
+
+	// Line12
+	auto line_pf_1 = SP::Ph1::PiLine::make("PiLine_13", Logger::Level::debug);
+	line_pf_1->setParameters(Three_bus_cfg.lineResistance, Three_bus_cfg.lineInductance, Three_bus_cfg.lineCapacitance);
+	line_pf_1->setBaseVoltage(Three_bus_cfg.Vnom);
+	// Line13
+	auto line_pf_2 = SP::Ph1::PiLine::make("PiLine_31", Logger::Level::debug);
+	line_pf_2->setParameters(Three_bus_cfg.lineResistance, Three_bus_cfg.lineInductance, Three_bus_cfg.lineCapacitance);
+	line_pf_2->setBaseVoltage(Three_bus_cfg.Vnom);
+	// Line23
+	auto line_pf_3 = SP::Ph1::PiLine::make("PiLine_23", Logger::Level::debug);
+	line_pf_3->setParameters(Three_bus_cfg.lineResistance, 2 * Three_bus_cfg.lineInductance, 2 * Three_bus_cfg.lineCapacitance);
+	line_pf_3->setBaseVoltage(Three_bus_cfg.Vnom);
+
+	// Topology
+	gen1_PF->connect({n1_PF});
+	gen2_PF->connect({n2_PF});
+	load_PF->connect({n3_PF});
+	line_pf_1->connect({n1_PF, n2_PF});
+	line_pf_2->connect({n1_PF, n3_PF});
+	line_pf_3->connect({n2_PF, n3_PF});
+
+	auto systemPF = SystemTopology(50,
+								   SystemNodeList{n1_PF, n2_PF, n3_PF},
+								   SystemComponentList{gen1_PF, gen2_PF, line_pf_1, line_pf_2, line_pf_3, load_PF});
+
+	//***Adding the PMU into the system topology
+	 systemPF.mComponents.push_back(pmu_1);
+	 systemPF.mComponents.push_back(pmu_2);
+	 systemPF.mComponents.push_back(pmu_3);
+
+	// Logging
+	auto loggerPF = DataLogger::make(simNamePF);
+	loggerPF->logAttribute("v_bus1", n1_PF->attribute("v"));
+	loggerPF->logAttribute("v_bus2", n2_PF->attribute("v"));
+	loggerPF->logAttribute("v_bus3", n3_PF->attribute("v"));
+
+	//***Logging the data from PMU
+	loggerPF->logAttribute("pmu_input_1", pmu_1->attribute("input"));
+	loggerPF->logAttribute("pmu_output_1", pmu_1->attribute("output"));
+
+	loggerPF->logAttribute("pmu_input_2", pmu_2->attribute("input"));
+	loggerPF->logAttribute("pmu_output_2", pmu_2->attribute("output"));
+	loggerPF->logAttribute("pmu_input_3", pmu_3->attribute("input"));
+	loggerPF->logAttribute("pmu_output_3", pmu_3->attribute("output"));
+
+	// Simulation
+	Simulation simPF(simNamePF, Logger::Level::debug);
+	simPF.setSystem(systemPF);
+	simPF.setTimeStep(timeStepPF);
+	simPF.setFinalTime(finalTimePF);
+	simPF.setDomain(Domain::SP);
+	simPF.setSolverType(Solver::Type::NRP);
+	simPF.setSolverAndComponentBehaviour(Solver::Behaviour::Initialization);
+	simPF.doInitFromNodesAndTerminals(false);
+	simPF.addLogger(loggerPF);
+
+	simPF.run();
+}
+
+int main(int argc, char *argv[])
+{
+
+	// Simultion parameters
+	String simName = "Three_bus_sim";
+	Real finalTime = 50;
+	Real timeStep = 0.001;
+	Real cmdInertia_G1 = 1.0;
+	Real cmdInertia_G2 = 1.0;
+	Real cmdDamping_G1 = 1.0;
+	Real cmdDamping_G2 = 1.0;
+
+	CommandLineArgs args(argc, argv);
+	if (argc > 1)
+	{
+		timeStep = args.timeStep;
+		finalTime = args.duration;
+		if (args.name != "dpsim")
+			simName = args.name;
+		if (args.options.find("SCALEINERTIA_G1") != args.options.end())
+			cmdInertia_G1 = args.getOptionReal("SCALEINERTIA_G1");
+		if (args.options.find("SCALEINERTIA_G2") != args.options.end())
+			cmdInertia_G2 = args.getOptionReal("SCALEINERTIA_G2");
+		if (args.options.find("SCALEDAMPING_G1") != args.options.end())
+			cmdDamping_G1 = args.getOptionReal("SCALEDAMPING_G1");
+		if (args.options.find("SCALEDAMPING_G2") != args.options.end())
+			cmdDamping_G2 = args.getOptionReal("SCALEDAMPING_G2");
+	}
+
+	Three_bus_sim(simName, timeStep, finalTime, cmdInertia_G1, cmdInertia_G2, cmdDamping_G1, cmdDamping_G2);
+}

dpsim/examples/cxx/Examples.h

@@ -472,6 +472,52 @@ namespace ThreeBus {
     };
 }
 
+namespace Three_bus_cfg{
+    struct ScenarioConfig {
+
+        //-----------Network-----------//
+        Real Vnom = 6.9e3;
+        Real nomFreq = 50;
+        Real nomOmega= nomFreq* 2*PI;
+
+        //-----------Generator 1 (bus1)-----------//
+        Real nomPower_G1 = 100e6;
+        Real nomPhPhVoltRMS_G1 = 6.9e3;
+        Real nomFreq_G1 = 50;
+
+        // Initialization parameters
+        Real initActivePower_G1 = 0.3386e6;
+        Real initReactivePower_G1= 0.3222e6;
+        Real setPointVoltage_G1=nomPhPhVoltRMS_G1;
+
+        //-----------Generator 2 (bus2)-----------//
+        Real nomPower_G2 = 100e6;
+        Real nomPhPhVoltRMS_G2 = 6.9e3;
+        Real nomFreq_G2 = 50;
+
+        // Initialization parameters
+        Real initActivePower_G2 = 0.2222e6;
+        Real initReactivePower_G2=(-1)*0.0285e6;
+        Real setPointVoltage_G2=nomPhPhVoltRMS_G2;
+
+        //-----------Transformers-----------//
+        Real t1_ratio=1;
+        Real t2_ratio=1;
+
+        //-----------Load (bus3)-----------
+        Real activePower_L= 0.5556e6;
+        Real reactivePower_L= 0.2778e6;
+
+        // -----------Transmission Lines-----------//
+
+        Real lineResistance = 0.025;
+        Real lineInductance = 0.075/(2*PI*50);
+        Real lineCapacitance = 0;
+
+    };
+
+}
+
 namespace SGIB {
 
     struct ScenarioConfig {

dpsim/src/PFSolver.cpp

@@ -65,10 +65,12 @@ void PFSolver::assignMatrixNodeIndices() {
 }
 
 void PFSolver::initializeComponents(){
-    for (auto comp : mSystem.mComponents) {
-        std::dynamic_pointer_cast<SimPowerComp<Complex>>(comp)->updateMatrixNodeIndices();
+   for (auto comp : mSystem.mComponents) {
+		auto sComp = std::dynamic_pointer_cast<SimSignalComp>(comp);
+		if (!sComp)	{
+		std::dynamic_pointer_cast<SimPowerComp<Complex>>(comp)->updateMatrixNodeIndices();
+		}
     }
-
 	SPDLOG_LOGGER_INFO(mSLog, "-- Initialize components from terminals or nodes of topology");
 	for (auto comp : mSystem.mComponents) {
 		auto pComp = std::dynamic_pointer_cast<SimPowerComp<Complex>>(comp);
@@ -404,5 +406,18 @@ void PFSolver::SolveTask::execute(Real time, Int timeStepCount) {
 }
 
 Task::List PFSolver::getTasks() {
-	return Task::List{std::make_shared<SolveTask>(*this)};
+	Task::List l;
+
+	for (auto comp : mSystem.mComponents) {
+	auto sComp = std::dynamic_pointer_cast<SimSignalComp>(comp);
+		if (sComp)	{
+		for (auto task : sComp->getTasks()) {
+			l.push_back(task);
+		}
+	}
+	}
+
+	l.push_back(std::make_shared<SolveTask>(*this));
+
+	return l;
 }