Commit 254bdb0e authored by Jan Dinkelbach's avatar Jan Dinkelbach Committed by Markus Mirz
Browse files

add emt3ph vs dp1ph comparison examples and notebook for vs, resistor and inductor

parent 68926ca6
......@@ -20,6 +20,7 @@ set(CIRCUIT_SOURCES
Circuits/DP_Slack_PiLine_VSI_with_PF_Init.cpp
Circuits/EMT_Slack_PiLine_PQLoad_with_PF_Init.cpp
Circuits/EMT_Slack_PiLine_VSI_with_PF_Init.cpp
Circuits/EMT_DP_VS_RLC_Circuits.cpp
# Powerflow examples
Circuits/PF_Slack_PiLine_PQLoad.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;
using namespace CPS::CIM;
void voltageSourceResistorEMT3ph() {
Real timeStep = 0.00005;
Real finalTime = 1;
String simName = "EMT_VoltageSource_Resistor";
Logger::setLogDir("logs/"+simName);
// Nodes
auto n1 = SimNode<Real>::make("n1", PhaseType::ABC);
// Components
auto vs = EMT::Ph3::VoltageSource::make("vs1");
vs->setParameters(CPS::Math::polar(100000, -PI/2.), 50);
auto res = EMT::Ph3::Resistor::make("R1");
res->setParameters(Reader::singlePhaseParameterToThreePhase(100));
// Topology
vs->connect({ SimNode<Real>::GND, n1 });
res->connect({n1, SimNode<Real>::GND});
auto sys = SystemTopology(50,
SystemNodeList{n1},
SystemComponentList{vs, res});
// Logging
auto logger = DataLogger::make(simName);
logger->addAttribute("v1", n1->attribute("v"));
logger->addAttribute("iR", res->attribute("i_intf"));
// Simulation
Simulation sim(simName);
sim.setSystem(sys);
sim.setTimeStep(timeStep);
sim.setFinalTime(finalTime);
sim.setDomain(Domain::EMT);
sim.addLogger(logger);
sim.run();
}
void voltageSourceResistorDP1ph() {
Real timeStep = 0.00005;
Real finalTime = 1;
String simName = "DP_VoltageSource_Resistor";
Logger::setLogDir("logs/"+simName);
// Nodes
auto n1 = SimNode<Complex>::make("n1");
// Components
auto vs = DP::Ph1::VoltageSource::make("vs1");
vs->setParameters(CPS::Math::polar(100000, -PI/2.));
auto res = DP::Ph1::Resistor::make("R1");
res->setParameters(100);
// Topology
vs->connect({ SimNode<Complex>::GND, n1 });
res->connect({n1, SimNode<Complex>::GND});
auto sys = SystemTopology(50,
SystemNodeList{n1},
SystemComponentList{vs, res});
// Logging
auto logger = DataLogger::make(simName);
logger->addAttribute("v1", n1->attribute("v"));
logger->addAttribute("iR", res->attribute("i_intf"));
// Simulation
Simulation sim(simName);
sim.setSystem(sys);
sim.setTimeStep(timeStep);
sim.setFinalTime(finalTime);
sim.setDomain(Domain::DP);
sim.addLogger(logger);
sim.run();
}
void voltageSourceInductorEMT3ph() {
Real timeStep = 0.00005;
Real finalTime = 1;
String simName = "EMT_VoltageSource_Inductor";
Logger::setLogDir("logs/"+simName);
// Nodes
auto n1 = SimNode<Real>::make("n1", PhaseType::ABC);
auto n2 = SimNode<Real>::make("n2", PhaseType::ABC);
// Components
auto vs = EMT::Ph3::VoltageSource::make("vs1", Logger::Level::debug);
vs->setParameters(CPS::Math::polar(100000, 0), 50);
auto res = EMT::Ph3::Resistor::make("R1");
res->setParameters(Reader::singlePhaseParameterToThreePhase(5));
auto ind = EMT::Ph3::Inductor::make("L1", Logger::Level::debug);
ind->setParameters(Reader::singlePhaseParameterToThreePhase(0.5));
// Topology
vs->connect({ SimNode<Real>::GND, n1 });
res->connect({n1, n2});
ind->connect({n2, SimNode<Real>::GND});
auto sys = SystemTopology(50,
SystemNodeList{n1, n2},
SystemComponentList{vs, res, ind});
// Logging
auto logger = DataLogger::make(simName);
logger->addAttribute("v1", n1->attribute("v"));
logger->addAttribute("v2", n2->attribute("v"));
logger->addAttribute("iR", res->attribute("i_intf"));
logger->addAttribute("iL", ind->attribute("i_intf"));
// Simulation
Simulation sim(simName);
sim.setSystem(sys);
sim.setTimeStep(timeStep);
sim.setFinalTime(finalTime);
sim.setDomain(Domain::EMT);
sim.addLogger(logger);
sim.run();
}
void voltageSourceInductorDP1ph() {
Real timeStep = 0.00005;
Real finalTime = 1;
String simName = "DP_VoltageSource_Inductor";
Logger::setLogDir("logs/"+simName);
// Nodes
auto n1 = SimNode<Complex>::make("n1");
auto n2 = SimNode<Complex>::make("n2");
// Components
auto vs = DP::Ph1::VoltageSource::make("vs1", Logger::Level::debug);
vs->setParameters(CPS::Math::polar(100000, 0));
auto res = DP::Ph1::Resistor::make("R1");
res->setParameters(5);
auto ind = DP::Ph1::Inductor::make("L1", Logger::Level::debug);
ind->setParameters(0.5);
// Topology
vs->connect({ SimNode<Complex>::GND, n1 });
res->connect({n1, n2});
ind->connect({n2, SimNode<Complex>::GND});
auto sys = SystemTopology(50,
SystemNodeList{n1, n2},
SystemComponentList{vs, res, ind});
// Logging
auto logger = DataLogger::make(simName);
logger->addAttribute("v1", n1->attribute("v"));
logger->addAttribute("v2", n2->attribute("v"));
logger->addAttribute("iR", res->attribute("i_intf"));
logger->addAttribute("iL", ind->attribute("i_intf"));
// Simulation
Simulation sim(simName);
sim.setSystem(sys);
sim.setTimeStep(timeStep);
sim.setFinalTime(finalTime);
sim.setDomain(Domain::DP);
sim.addLogger(logger);
sim.run();
}
int main(int argc, char* argv[]) {
voltageSourceResistorEMT3ph();
voltageSourceResistorDP1ph();
voltageSourceInductorEMT3ph();
voltageSourceInductorDP1ph();
}
%% Cell type:code id: tags:
``` python
%%bash
TOP=${TOP:-$(git rev-parse --show-toplevel)}
PATH=${TOP}/build/Examples/Cxx
EMT_DP_VS_RLC_Circuits
```
%% Cell type:code id: tags:
``` python
from villas.dataprocessing.readtools import *
from villas.dataprocessing.timeseries import *
from villas.dataprocessing.timeseries import TimeSeries as ts
import matplotlib.pyplot as plt
import numpy as np
import re
%matplotlib widget
```
%% Cell type:markdown id: tags:
## Voltage Source + Resistor
%% Cell type:code id: tags:
``` python
model_name = 'VoltageSource_Resistor'
path_DP = 'logs/' + 'DP_' + model_name + '/'
dpsim_result_file_DP = path_DP + 'DP_' + model_name + '.csv'
ts_dpsim_DP = read_timeseries_csv(dpsim_result_file_DP)
ts_dpsim_DP_shifted = ts.frequency_shift_list(ts_dpsim_DP, 50)
path_EMT = 'logs/' + 'EMT_' + model_name + '/'
dpsim_result_file_EMT = path_EMT + 'EMT_' + model_name + '.csv'
ts_dpsim_EMT = read_timeseries_csv(dpsim_result_file_EMT)
```
%% Cell type:markdown id: tags:
### Plot voltages
%% Cell type:code id: tags:
``` python
plt.figure()
# EMT
var_names = ['v1_0']
for var_name in var_names:
plt.plot(ts_dpsim_EMT[var_name].time, np.sqrt(3/2)*ts_dpsim_EMT[var_name].values, label=var_name)
# DP
var_names = ['v1_shift']
for var_name in var_names:
plt.plot(ts_dpsim_DP_shifted[var_name].time, ts_dpsim_DP_shifted[var_name].values, label=var_name, linestyle=':')
plt.legend()
plt.show()
```
%% Cell type:markdown id: tags:
### Plot currents
%% Cell type:code id: tags:
``` python
plt.figure()
# EMT
var_names = ['iR_0']
for var_name in var_names:
plt.plot(ts_dpsim_EMT[var_name].time, np.sqrt(3/2)*ts_dpsim_EMT[var_name].values, label=var_name)
# DP
var_names = ['iR_shift']
for var_name in var_names:
plt.plot(ts_dpsim_DP_shifted[var_name].time, ts_dpsim_DP_shifted[var_name].values, label=var_name, linestyle=':')
plt.legend()
plt.show()
```
%% Cell type:markdown id: tags:
### Comparison DP vs. EMT
%% Cell type:code id: tags:
``` python
plt.figure()
for name in [('v1_0', 'v1_shift')]:
plt.plot(ts_dpsim_EMT[name[0]].time, np.sqrt(3/2)*ts_dpsim_EMT[name[0]].values - ts_dpsim_DP_shifted[name[1]].values, label=name[0]+' vs. '+name[1])
plt.legend()
```
%% Cell type:code id: tags:
``` python
plt.figure()
for name in [('iR_0', 'iR_shift')]:
plt.plot(ts_dpsim_EMT[name[0]].time, np.sqrt(3/2)*ts_dpsim_EMT[name[0]].values - ts_dpsim_DP_shifted[name[1]].values, label=name[0]+' vs. '+name[1])
plt.legend()
```
%% Cell type:markdown id: tags:
### Assertion DP vs. EMT
%% Cell type:code id: tags:
``` python
compare_errors_abs = []
compare_errors_rel = []
for name in [('v1_0', 'v1_shift'), ('iR_0', 'iR_shift')]:
compare_errors_abs.append(np.absolute(np.sqrt(3/2)*ts_dpsim_EMT[name[0]].values - ts_dpsim_DP_shifted[name[1]].values).max())
compare_errors_rel.append(np.absolute(np.sqrt(3/2)*ts_dpsim_EMT[name[0]].values - ts_dpsim_DP_shifted[name[1]].values).max()/ts_dpsim_DP_shifted[name[1]].values.max())
print(name[0]+' vs. '+name[1] + ' (abs): ' + str(compare_errors_abs[-1]))
print(name[0]+' vs. '+name[1] + ' (rel): ' + str(compare_errors_rel[-1]))
print('Max rel error: '+ '{:.2}'.format(np.max(compare_errors_rel)*100) +'%')
assert np.max(compare_errors_rel) < 1e-3
```
%% Cell type:markdown id: tags:
## Voltage Source + Resistor + Inductor
%% Cell type:code id: tags:
``` python
model_name = 'VoltageSource_Inductor'
path_DP = 'logs/' + 'DP_' + model_name + '/'
dpsim_result_file_DP = path_DP + 'DP_' + model_name + '.csv'
ts_dpsim_DP = read_timeseries_csv(dpsim_result_file_DP)
ts_dpsim_DP_shifted = ts.frequency_shift_list(ts_dpsim_DP, 50)
path_EMT = 'logs/' + 'EMT_' + model_name + '/'
dpsim_result_file_EMT = path_EMT + 'EMT_' + model_name + '.csv'
ts_dpsim_EMT = read_timeseries_csv(dpsim_result_file_EMT)
```
%% Cell type:markdown id: tags:
### Plot voltages
%% Cell type:code id: tags:
``` python
plt.figure()
# EMT
var_names = ['v1_0', 'v2_0']
for var_name in var_names:
plt.plot(ts_dpsim_EMT[var_name].time, np.sqrt(3/2)*ts_dpsim_EMT[var_name].values, label=var_name)
# DP
var_names = ['v1_shift', 'v2_shift']
for var_name in var_names:
plt.plot(ts_dpsim_DP_shifted[var_name].time, ts_dpsim_DP_shifted[var_name].values, label=var_name, linestyle=':')
plt.legend()
plt.show()
```
%% Cell type:markdown id: tags:
### Plot currents
%% Cell type:code id: tags:
``` python
plt.figure()
# EMT
var_names = ['iR_0'] # ['iR_0', 'iL_0']
for var_name in var_names:
plt.plot(ts_dpsim_EMT[var_name].time, np.sqrt(3/2)*ts_dpsim_EMT[var_name].values, label=var_name)
# DP
var_names = ['iR_shift'] # ['iR_shift', 'iL_shift']
for var_name in var_names:
plt.plot(ts_dpsim_DP_shifted[var_name].time, ts_dpsim_DP_shifted[var_name].values, label=var_name, linestyle=':')
plt.legend()
plt.show()
```
%% Cell type:markdown id: tags:
### Comparison DP vs. EMT
%% Cell type:code id: tags:
``` python
plt.figure()
for name in [('v1_0', 'v1_shift'), ('v2_0', 'v2_shift')]:
plt.plot(ts_dpsim_EMT[name[0]].time, np.sqrt(3/2)*ts_dpsim_EMT[name[0]].values - ts_dpsim_DP_shifted[name[1]].values, label=name[0]+' vs. '+name[1])
plt.legend()
```
%% Cell type:code id: tags:
``` python
plt.figure()
for name in [('iR_0', 'iR_shift'), ('iL_0', 'iL_shift')]:
plt.plot(ts_dpsim_EMT[name[0]].time, np.sqrt(3/2)*ts_dpsim_EMT[name[0]].values - ts_dpsim_DP_shifted[name[1]].values, label=name[0]+' vs. '+name[1])
plt.legend()
```
%% Cell type:markdown id: tags:
### Assertion DP vs. EMT
%% Cell type:code id: tags:
``` python
compare_errors_abs = []
compare_errors_rel = []
for name in [('v1_0', 'v1_shift'), ('v2_0', 'v2_shift'), ('iR_0', 'iR_shift'), ('iL_0', 'iL_shift')]:
compare_errors_abs.append(np.absolute(np.sqrt(3/2)*ts_dpsim_EMT[name[0]].values - ts_dpsim_DP_shifted[name[1]].values).max())
compare_errors_rel.append(np.absolute(np.sqrt(3/2)*ts_dpsim_EMT[name[0]].values - ts_dpsim_DP_shifted[name[1]].values).max()/ts_dpsim_DP_shifted[name[1]].values.max())
print(name[0]+' vs. '+name[1] + ' (abs): ' + str(compare_errors_abs[-1]))
print(name[0]+' vs. '+name[1] + ' (rel): ' + str(compare_errors_rel[-1]))
print('Max rel error: '+ '{:.2}'.format(np.max(compare_errors_rel)*100) +'%')
assert np.max(compare_errors_rel) < 1e-3
```
%% Cell type:code id: tags:
``` python
```
......@@ -73,9 +73,15 @@ void EMT::Ph3::Inductor::mnaInitialize(Real omega, Real timeStep, Attribute<Matr
mMnaTasks.push_back(std::make_shared<MnaPostStep>(*this, leftVector));
mRightVector = Matrix::Zero(leftVector->get().rows(), 1);
mSLog->debug(
"\nEquivalent Current (mnaInitialize): {:s}",
Logger::matrixToString(mEquivCurrent));
mSLog->info(
"\n--- MNA initialization ---"
"\nInitial voltage {:s}"
"\nInitial current {:s}"
"\nEquiv. current {:s}"
"\n--- MNA initialization finished ---",
Logger::matrixToString(mIntfVoltage),
Logger::matrixToString(mIntfCurrent),
Logger::matrixToString(mEquivCurrent));
mSLog->flush();
}
......@@ -187,12 +193,16 @@ void EMT::Ph3::Inductor::mnaUpdateVoltage(const Matrix& leftVector) {
mIntfVoltage(1, 0) = mIntfVoltage(1, 0) - Math::realFromVectorElement(leftVector, matrixNodeIndex(0, 1));
mIntfVoltage(2, 0) = mIntfVoltage(2, 0) - Math::realFromVectorElement(leftVector, matrixNodeIndex(0, 2));
}
mSLog->debug(
"\nUpdate Voltage: {:s}",
Logger::matrixToString(mIntfVoltage)
);
}
void EMT::Ph3::Inductor::mnaUpdateCurrent(const Matrix& leftVector) {
mIntfCurrent = mEquivCond * mIntfVoltage + mEquivCurrent;
mSLog->debug(
"\nCurrent: {:s}",
"\nUpdate Current: {:s}",
Logger::matrixToString(mIntfCurrent)
);
mSLog->flush();
......
Markdown is supported
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