Commit 39acad66 authored by Jan Dinkelbach's avatar Jan Dinkelbach Committed by Markus Mirz
Browse files

update emt3ph network injection and add validation example

parent 69c26264
......@@ -17,27 +17,28 @@
#include <DPsim.h>
using namespace DPsim;
using namespace CPS::DP;
using namespace CPS;
using namespace CPS::CIM;
void simElements() {
void simElementsDP1ph() {
Real timeStep = 0.00005;
Real finalTime = 1;
String simName = "DP_Slack_Elements";
Logger::setLogDir("logs/"+simName);
// Nodes
auto n1 = SimNode::make("n1");
auto n1 = SimNode<Complex>::SimNode::make("n1");
// Components
auto vs = Ph1::VoltageSource::make("v_1");
auto vs = DP::Ph1::VoltageSource::make("v_1");
vs->setParameters(CPS::Math::polar(100000, 0));
auto load = Ph1::Resistor::make("R_load");
auto load = DP::Ph1::Resistor::make("R_load");
load->setParameters(10000);
// Topology
vs->connect({ SimNode::GND, n1 });
load->connect({ n1, SimNode::GND });
vs->connect({ SimNode<Complex>::GND, n1 });
load->connect({ n1, SimNode<Complex>::GND });
auto sys = SystemTopology(50,
SystemNodeList{n1},
......@@ -57,25 +58,25 @@ void simElements() {
sim.run();
}
void simComponent() {
void simComponentDP1ph() {
Real timeStep = 0.00005;
Real finalTime = 1;
String simName = "DP_Slack_Component";
Logger::setLogDir("logs/"+simName);
// Nodes
auto n1 = SimNode::make("n1");
auto n1 = SimNode<Complex>::make("n1");
// Components
auto sl = Ph1::NetworkInjection::make("v_1");
auto sl = DP::Ph1::NetworkInjection::make("v_1");
sl->setParameters(CPS::Math::polar(100000, 0));
auto load = Ph1::Resistor::make("R_load");
auto load = DP::Ph1::Resistor::make("R_load");
load->setParameters(10000);
// Topology
sl->connect({ n1 });
load->connect({ n1, SimNode::GND });
load->connect({ n1, SimNode<Complex>::GND });
auto sys = SystemTopology(50,
SystemNodeList{n1},
......@@ -95,7 +96,88 @@ void simComponent() {
sim.run();
}
void simElementsEMT3ph() {
Real timeStep = 0.00005;
Real finalTime = 1;
String simName = "EMT_Slack_Elements";
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, 0), 50);
auto load = EMT::Ph3::Resistor::make("Rload");
load->setParameters(Reader::singlePhaseParameterToThreePhase(10000));
// Topology
vs->connect({ SimNode<Real>::GND, n1 });
load->connect({n1, SimNode<Real>::GND});
auto sys = SystemTopology(50,
SystemNodeList{n1},
SystemComponentList{vs, load});
// Logging
auto logger = DataLogger::make(simName);
logger->addAttribute("v1", n1->attribute("v"));
logger->addAttribute("i1", vs->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 simComponentEMT3ph() {
Real timeStep = 0.00005;
Real finalTime = 1;
String simName = "EMT_Slack_Component";
Logger::setLogDir("logs/"+simName);
// Nodes
auto n1 = SimNode<Real>::make("n1", PhaseType::ABC);
// Components
auto vs = EMT::Ph3::NetworkInjection::make("vs1", Logger::Level::debug);
vs->setParameters(Reader::singlePhaseVariableToThreePhase(CPS::Math::polar(100000, 0)), 50);
auto load = EMT::Ph3::Resistor::make("Rload");
load->setParameters(Reader::singlePhaseParameterToThreePhase(10000));
// Topology
vs->connect({ n1 });
load->connect({n1, SimNode<Real>::GND});
auto sys = SystemTopology(50,
SystemNodeList{n1},
SystemComponentList{vs, load});
// Logging
auto logger = DataLogger::make(simName);
logger->addAttribute("v1", n1->attribute("v"));
logger->addAttribute("i1", vs->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();
}
int main(int argc, char* argv[]) {
simElements();
simComponent();
simElementsDP1ph();
simComponentDP1ph();
simElementsEMT3ph();
simComponentEMT3ph();
}
\ No newline at end of file
......@@ -114,5 +114,155 @@
for name in ['v1_shift', 'i1_shift']:
errors_dp_shifted.append(np.absolute(slack_elements_dp_shifted[name].values - slack_component_dp_shifted[name].values).max())
print(name + ': ' + str(errors_dp_shifted[-1]))
assert np.max(errors_dp_shifted) < 1e-3
```
%% Cell type:markdown id: tags:
## EMT Slack with voltage source element
%% Cell type:code id: tags:
``` python
work_dir = 'logs/EMT_Slack_Elements/'
log_name = 'EMT_Slack_Elements'
print(work_dir + log_name + '.csv')
slack_elements_emt = rt.read_timeseries_dpsim(work_dir + log_name + '.csv')
```
%% Cell type:code id: tags:
``` python
plt.figure()
plt.plot(slack_elements_emt['v1_0'].time, slack_elements_emt['v1_0'].values, label='v1_0')
plt.plot(slack_elements_emt['v1_1'].time, slack_elements_emt['v1_1'].values, label='v1_1')
plt.plot(slack_elements_emt['v1_2'].time, slack_elements_emt['v1_2'].values, label='v1_2')
plt.legend()
```
%% Cell type:code id: tags:
``` python
plt.figure()
plt.plot(slack_elements_emt['i1_0'].time, slack_elements_emt['i1_0'].values, label='i1_0')
plt.plot(slack_elements_emt['i1_1'].time, slack_elements_emt['i1_1'].values, label='i1_1')
plt.plot(slack_elements_emt['i1_2'].time, slack_elements_emt['i1_2'].values, label='i1_2')
plt.legend()
```
%% Cell type:markdown id: tags:
## EMT Slack composite model
%% Cell type:code id: tags:
``` python
work_dir = 'logs/EMT_Slack_Component/'
log_name = 'EMT_Slack_Component'
print(work_dir + log_name + '.csv')
slack_component_emt = rt.read_timeseries_dpsim(work_dir + log_name + '.csv')
```
%% Cell type:code id: tags:
``` python
plt.figure()
plt.plot(slack_component_emt['v1_0'].time, slack_component_emt['v1_0'].values, label='v1_0')
plt.plot(slack_component_emt['v1_1'].time, slack_component_emt['v1_1'].values, label='v1_1')
plt.plot(slack_component_emt['v1_2'].time, slack_component_emt['v1_2'].values, label='v1_2')
plt.legend()
```
%% Cell type:code id: tags:
``` python
plt.figure()
plt.plot(slack_component_emt['i1_0'].time, slack_component_emt['i1_0'].values, label='i1_0')
plt.plot(slack_component_emt['i1_1'].time, slack_component_emt['i1_1'].values, label='i1_1')
plt.plot(slack_component_emt['i1_2'].time, slack_component_emt['i1_2'].values, label='i1_2')
plt.legend()
```
%% Cell type:markdown id: tags:
## Error for EMT Slack
%% Cell type:code id: tags:
``` python
plt.figure()
plt.plot(slack_elements_emt['v1_0'].time, slack_elements_emt['v1_0'].values - slack_component_emt['v1_0'].values, label='v1_0_error')
plt.plot(slack_elements_emt['v1_1'].time, slack_elements_emt['v1_1'].values - slack_component_emt['v1_1'].values, label='v1_1_error')
plt.plot(slack_elements_emt['v1_2'].time, slack_elements_emt['v1_2'].values - slack_component_emt['v1_2'].values, label='v1_2_error')
plt.legend()
```
%% Cell type:code id: tags:
``` python
plt.figure()
plt.plot(slack_elements_emt['i1_0'].time, slack_elements_emt['i1_0'].values - slack_component_emt['i1_0'].values, label='i1_0_error')
plt.plot(slack_elements_emt['i1_1'].time, slack_elements_emt['i1_1'].values - slack_component_emt['i1_1'].values, label='i1_1_error')
plt.plot(slack_elements_emt['i1_2'].time, slack_elements_emt['i1_2'].values - slack_component_emt['i1_2'].values, label='i1_2_error')
plt.legend()
```
%% Cell type:markdown id: tags:
## Assertion for EMT Slack
%% Cell type:code id: tags:
``` python
errors_emt = []
for name in ['v1_0', 'v1_1', 'v1_2', 'i1_0', 'i1_1', 'i1_2']:
errors_emt.append(np.absolute(slack_elements_emt[name].values - slack_component_emt[name].values).max())
print(name + ': ' + str(errors_emt[-1]))
assert np.max(errors_emt) < 1e-3
```
%% 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(slack_component_emt[name[0]].time, np.sqrt(3/2)*slack_component_emt[name[0]].values - slack_component_dp_shifted[name[1]].values, label=name[0]+' vs. '+name[1])
plt.legend()
```
%% Cell type:code id: tags:
``` python
plt.figure()
for name in [('i1_0', 'i1_shift')]:
plt.plot(slack_component_emt[name[0]].time, np.sqrt(3/2)*slack_component_emt[name[0]].values - slack_component_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'), ('i1_0', 'i1_shift')]:
compare_errors_abs.append(np.absolute(np.sqrt(3/2)*slack_component_emt[name[0]].values - slack_component_dp_shifted[name[1]].values).max())
compare_errors_rel.append(np.absolute(np.sqrt(3/2)*slack_component_emt[name[0]].values - slack_component_dp_shifted[name[1]].values).max()/slack_component_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
```
......
......@@ -178,6 +178,8 @@ namespace CIM {
///
void initDynamicSystemTopologyWithPowerflow(SystemTopology& systemPF, SystemTopology& system);
/// To convert single phase complex variables (voltages, currents) to symmetrical three phase ones
static MatrixComp singlePhaseVariableToThreePhase(Complex var_1ph);
///
static Matrix singlePhaseParameterToThreePhase(Real parameter);
///
......
......@@ -86,6 +86,15 @@ TopologicalPowerComp::Ptr Reader::mapComponent(BaseClass* obj) {
return nullptr;
}
MatrixComp Reader::singlePhaseVariableToThreePhase(Complex var_1ph) {
MatrixComp var_3ph = MatrixComp::Zero(3, 1);
var_3ph <<
var_1ph,
var_1ph * SHIFT_TO_PHASE_B,
var_1ph * SHIFT_TO_PHASE_C;
return var_3ph;
}
Matrix Reader::singlePhaseParameterToThreePhase(Real parameter) {
Matrix param_3ph = Matrix::Zero(3, 3);
param_3ph <<
......
......@@ -44,20 +44,26 @@ void EMT::Ph3::NetworkInjection::initializeFromNodesAndTerminals(Real frequency)
mSrcFreq = attribute<Real>("f_src");
mSrcFreq->set(frequency);
MatrixComp vInitABC = MatrixComp::Zero(3, 1);
vInitABC(0, 0) = RMS3PH_TO_PEAK1PH * initialSingleVoltage(0);
vInitABC(1, 0) = vInitABC(0, 0) * SHIFT_TO_PHASE_B;
vInitABC(2, 0) = vInitABC(0, 0) * SHIFT_TO_PHASE_C;
mVoltageRef->set(vInitABC);
mSLog->info(
"\n--- Initialization from node voltages ---"
"\nReference voltage: {:s}"
"\nTerminal 0 voltage: {:s}"
"\n--- Initialization from node voltages ---",
Logger::matrixCompToString(mVoltageRef->get()),
Logger::phasorToString(RMS3PH_TO_PEAK1PH * initialSingleVoltage(0)));
mSLog->info("\n--- Initialization from node voltages ---");
// TODO: this approach currently does not work, if voltage ref set from outside without using setParameters,
// since mParametersSet remains false then
if (!mParametersSet) {
MatrixComp vInitABC = MatrixComp::Zero(3, 1);
vInitABC(0, 0) = initialSingleVoltage(0);
vInitABC(1, 0) = initialSingleVoltage(0) * SHIFT_TO_PHASE_B;
vInitABC(2, 0) = initialSingleVoltage(0) * SHIFT_TO_PHASE_C;
mVoltageRef->set(vInitABC);
mSLog->info("\nReference voltage: {:s}"
"\nTerminal 0 voltage: {:s}",
Logger::matrixCompToString(mVoltageRef->get()),
Logger::phasorToString(initialSingleVoltage(0)));
} else {
mSLog->info("\nInitialization from node voltages omitted (parameter already set)."
"\nReference voltage: {:s}",
Logger::matrixCompToString(mVoltageRef->get()));
}
mSLog->info("\n--- Initialization from node voltages ---");
mSLog->flush();
}
......@@ -111,15 +117,15 @@ void EMT::Ph3::NetworkInjection::mnaApplyRightSideVectorStamp(Matrix& rightVecto
void EMT::Ph3::NetworkInjection::updateVoltage(Real time) {
if (mSrcFreq->get() < 0) {
mIntfVoltage = mVoltageRef->get().real();
mIntfVoltage = RMS3PH_TO_PEAK1PH * mVoltageRef->get().real();
}
else {
mIntfVoltage(0, 0) =
Math::abs(mVoltageRef->get()(0, 0)) * cos(time * 2. * PI * mSrcFreq->get() + Math::phase(mVoltageRef->get())(0, 0));
RMS3PH_TO_PEAK1PH * Math::abs(mVoltageRef->get()(0, 0)) * cos(time * 2. * PI * mSrcFreq->get() + Math::phase(mVoltageRef->get())(0, 0));
mIntfVoltage(1, 0) =
Math::abs(mVoltageRef->get()(1, 0)) * cos(time * 2. * PI * mSrcFreq->get() + Math::phase(mVoltageRef->get())(1, 0));
RMS3PH_TO_PEAK1PH * Math::abs(mVoltageRef->get()(1, 0)) * cos(time * 2. * PI * mSrcFreq->get() + Math::phase(mVoltageRef->get())(1, 0));
mIntfVoltage(2, 0) =
Math::abs(mVoltageRef->get()(2, 0)) * cos(time * 2. * PI * mSrcFreq->get() + Math::phase(mVoltageRef->get())(2, 0));
RMS3PH_TO_PEAK1PH * Math::abs(mVoltageRef->get()(2, 0)) * cos(time * 2. * PI * mSrcFreq->get() + Math::phase(mVoltageRef->get())(2, 0));
}
}
......
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