DP_Ph1_VoltageSource.cpp 9.15 KB
Newer Older
1
2
/* Copyright 2017-2020 Institute for Automation of Complex Power Systems,
 *                     EONERC, RWTH Aachen University
Markus Mirz's avatar
Markus Mirz committed
3
 *
4
5
6
 * 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/.
Markus Mirz's avatar
Markus Mirz committed
7
8
9
10
11
12
13
 *********************************************************************************/

#include <cps/DP/DP_Ph1_VoltageSource.h>

using namespace CPS;

DP::Ph1::VoltageSource::VoltageSource(String uid, String name, Logger::Level logLevel)
14
	: SimPowerComp<Complex>(uid, name, logLevel) {
Markus Mirz's avatar
Markus Mirz committed
15
16
17
18
19
20
21
22
23
	setVirtualNodeNumber(1);
	setTerminalNumber(2);
	mIntfVoltage = MatrixComp::Zero(1,1);
	mIntfCurrent = MatrixComp::Zero(1,1);

	addAttribute<Complex>("V_ref", Flags::read | Flags::write);
	addAttribute<Real>("f_src", Flags::read | Flags::write);
}

Markus Mirz's avatar
Markus Mirz committed
24
SimPowerComp<Complex>::Ptr DP::Ph1::VoltageSource::clone(String name) {
Markus Mirz's avatar
Markus Mirz committed
25
26
27
28
29
30
	auto copy = VoltageSource::make(name, mLogLevel);
	copy->setParameters(attribute<Complex>("V_ref")->get());
	return copy;
}

void DP::Ph1::VoltageSource::initialize(Matrix frequencies) {
Markus Mirz's avatar
Markus Mirz committed
31
	SimPowerComp<Complex>::initialize(frequencies);
Markus Mirz's avatar
Markus Mirz committed
32
33
34
35
36
}

void DP::Ph1::VoltageSource::setParameters(Complex voltageRef, Real srcFreq) {
	attribute<Complex>("V_ref")->set(voltageRef);
	attribute<Real>("f_src")->set(srcFreq);
37

38
	mParametersSet = true;
Markus Mirz's avatar
Markus Mirz committed
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
}

void DP::Ph1::VoltageSource::initializeFromPowerflow(Real frequency) {
	mVoltageRef = attribute<Complex>("V_ref");
	mSrcFreq = attribute<Real>("f_src");
	if (mVoltageRef->get() == Complex(0, 0))
		mVoltageRef->set(initialSingleVoltage(1) - initialSingleVoltage(0));

	mSLog->info(
		"\n--- Initialization from node voltages ---"
		"\nVoltage across: {:e}<{:e}"
		"\nTerminal 0 voltage: {:e}<{:e}"
		"\nTerminal 1 voltage: {:e}<{:e}"
		"\n--- Initialization from node voltages ---",
		std::abs(mVoltageRef->get()), std::arg(mVoltageRef->get()),
		std::abs(initialSingleVoltage(0)), std::arg(initialSingleVoltage(0)),
		std::abs(initialSingleVoltage(1)), std::arg(initialSingleVoltage(1)));
}

// #### MNA functions ####

60
61
62
63
64
65
66
67
68
69
70
void DP::Ph1::VoltageSource::mnaAddPreStepDependencies(AttributeBase::List &prevStepDependencies, AttributeBase::List &attributeDependencies, AttributeBase::List &modifiedAttributes) {
	attributeDependencies.push_back(attribute("V_ref"));
	modifiedAttributes.push_back(attribute("right_vector"));
	modifiedAttributes.push_back(attribute("v_intf"));
}

void DP::Ph1::VoltageSource::mnaAddPostStepDependencies(AttributeBase::List &prevStepDependencies, AttributeBase::List &attributeDependencies, AttributeBase::List &modifiedAttributes, Attribute<Matrix>::Ptr &leftVector) {
	attributeDependencies.push_back(leftVector);
	modifiedAttributes.push_back(attribute("i_intf"));
};

Markus Mirz's avatar
Markus Mirz committed
71
72
void DP::Ph1::VoltageSource::mnaInitialize(Real omega, Real timeStep, Attribute<Matrix>::Ptr leftVector) {
	MNAInterface::mnaInitialize(omega, timeStep);
Markus Mirz's avatar
Markus Mirz committed
73
	updateMatrixNodeIndices();
Markus Mirz's avatar
Markus Mirz committed
74
75
76
77
78
79
80
81
82

	mIntfVoltage(0,0) = mVoltageRef->get();
	mMnaTasks.push_back(std::make_shared<MnaPreStep>(*this));
	mMnaTasks.push_back(std::make_shared<MnaPostStep>(*this, leftVector));
	mRightVector = Matrix::Zero(leftVector->get().rows(), 1);
}

void DP::Ph1::VoltageSource::mnaInitializeHarm(Real omega, Real timeStep, std::vector<Attribute<Matrix>::Ptr> leftVectors) {
	MNAInterface::mnaInitialize(omega, timeStep);
Markus Mirz's avatar
Markus Mirz committed
83
	updateMatrixNodeIndices();
Markus Mirz's avatar
Markus Mirz committed
84
85
86
87
88
89
90
91
92
93
94

	mIntfVoltage(0,0) = mVoltageRef->get();

	mMnaTasks.push_back(std::make_shared<MnaPreStepHarm>(*this));
	mMnaTasks.push_back(std::make_shared<MnaPostStepHarm>(*this, leftVectors));
	mRightVector = Matrix::Zero(leftVectors[0]->get().rows(), mNumFreqs);
}

void DP::Ph1::VoltageSource::mnaApplySystemMatrixStamp(Matrix& systemMatrix) {
	for (UInt freq = 0; freq < mNumFreqs; freq++) {
		if (terminalNotGrounded(0)) {
Markus Mirz's avatar
Markus Mirz committed
95
96
			Math::setMatrixElement(systemMatrix, mVirtualNodes[0]->matrixNodeIndex(), matrixNodeIndex(0), Complex(-1, 0), mNumFreqs, freq);
			Math::setMatrixElement(systemMatrix, matrixNodeIndex(0), mVirtualNodes[0]->matrixNodeIndex(), Complex(-1, 0), mNumFreqs, freq);
Markus Mirz's avatar
Markus Mirz committed
97
98
		}
		if (terminalNotGrounded(1)) {
Markus Mirz's avatar
Markus Mirz committed
99
100
			Math::setMatrixElement(systemMatrix, mVirtualNodes[0]->matrixNodeIndex(), matrixNodeIndex(1), Complex(1, 0), mNumFreqs, freq);
			Math::setMatrixElement(systemMatrix, matrixNodeIndex(1), mVirtualNodes[0]->matrixNodeIndex(), Complex(1, 0), mNumFreqs, freq);
Markus Mirz's avatar
Markus Mirz committed
101
102
103
104
		}

		mSLog->info("-- Stamp frequency {:d} ---", freq);
		if (terminalNotGrounded(0)) {
Markus Mirz's avatar
Markus Mirz committed
105
106
			mSLog->info("Add {:f} to system at ({:d},{:d})", -1., matrixNodeIndex(0), mVirtualNodes[0]->matrixNodeIndex());
			mSLog->info("Add {:f} to system at ({:d},{:d})", -1., mVirtualNodes[0]->matrixNodeIndex(), matrixNodeIndex(0));
Markus Mirz's avatar
Markus Mirz committed
107
108
		}
		if (terminalNotGrounded(1)) {
Markus Mirz's avatar
Markus Mirz committed
109
110
			mSLog->info("Add {:f} to system at ({:d},{:d})", 1., mVirtualNodes[0]->matrixNodeIndex(), matrixNodeIndex(1));
			mSLog->info("Add {:f} to system at ({:d},{:d})", 1., matrixNodeIndex(1), mVirtualNodes[0]->matrixNodeIndex());
Markus Mirz's avatar
Markus Mirz committed
111
112
113
114
115
116
		}
	}
}

void DP::Ph1::VoltageSource::mnaApplySystemMatrixStampHarm(Matrix& systemMatrix, Int freqIdx) {
	if (terminalNotGrounded(0)) {
Markus Mirz's avatar
Markus Mirz committed
117
118
		Math::setMatrixElement(systemMatrix, mVirtualNodes[0]->matrixNodeIndex(), matrixNodeIndex(0), Complex(-1, 0));
		Math::setMatrixElement(systemMatrix, matrixNodeIndex(0), mVirtualNodes[0]->matrixNodeIndex(), Complex(-1, 0));
Markus Mirz's avatar
Markus Mirz committed
119
120
	}
	if (terminalNotGrounded(1)) {
Markus Mirz's avatar
Markus Mirz committed
121
122
		Math::setMatrixElement(systemMatrix, mVirtualNodes[0]->matrixNodeIndex(), matrixNodeIndex(1), Complex(1, 0));
		Math::setMatrixElement(systemMatrix, matrixNodeIndex(1), mVirtualNodes[0]->matrixNodeIndex(), Complex(1, 0));
Markus Mirz's avatar
Markus Mirz committed
123
124
125
126
	}

	mSLog->info("-- Stamp frequency {:d} ---", freqIdx);
	if (terminalNotGrounded(0)) {
Markus Mirz's avatar
Markus Mirz committed
127
128
		mSLog->info("Add {:f} to system at ({:d},{:d})", -1., matrixNodeIndex(0), mVirtualNodes[0]->matrixNodeIndex());
		mSLog->info("Add {:f} to system at ({:d},{:d})", -1., mVirtualNodes[0]->matrixNodeIndex(), matrixNodeIndex(0));
Markus Mirz's avatar
Markus Mirz committed
129
130
	}
	if (terminalNotGrounded(1)) {
Markus Mirz's avatar
Markus Mirz committed
131
132
		mSLog->info("Add {:f} to system at ({:d},{:d})", 1., mVirtualNodes[0]->matrixNodeIndex(), matrixNodeIndex(1));
		mSLog->info("Add {:f} to system at ({:d},{:d})", 1., matrixNodeIndex(1), mVirtualNodes[0]->matrixNodeIndex());
Markus Mirz's avatar
Markus Mirz committed
133
134
135
136
137
	}
}

void DP::Ph1::VoltageSource::mnaApplyRightSideVectorStamp(Matrix& rightVector) {
	// TODO: Is this correct with two nodes not gnd?
Markus Mirz's avatar
Markus Mirz committed
138
	Math::setVectorElement(rightVector, mVirtualNodes[0]->matrixNodeIndex(), mIntfVoltage(0,0), mNumFreqs);
Markus Mirz's avatar
Markus Mirz committed
139
	SPDLOG_LOGGER_DEBUG(mSLog, "Add {:s} to source vector at {:d}",
Markus Mirz's avatar
Markus Mirz committed
140
		Logger::complexToString(mIntfVoltage(0,0)), mVirtualNodes[0]->matrixNodeIndex());
Markus Mirz's avatar
Markus Mirz committed
141
142
143
144
145
}

void DP::Ph1::VoltageSource::mnaApplyRightSideVectorStampHarm(Matrix& rightVector) {
	for (UInt freq = 0; freq < mNumFreqs; freq++) {
		// TODO: Is this correct with two nodes not gnd?
Markus Mirz's avatar
Markus Mirz committed
146
		Math::setVectorElement(rightVector, mVirtualNodes[0]->matrixNodeIndex(), mIntfVoltage(0,freq), 1, 0, freq);
Markus Mirz's avatar
Markus Mirz committed
147
		SPDLOG_LOGGER_DEBUG(mSLog, "Add {:s} to source vector at {:d}",
Markus Mirz's avatar
Markus Mirz committed
148
			Logger::complexToString(mIntfVoltage(0,freq)), mVirtualNodes[0]->matrixNodeIndex());
Markus Mirz's avatar
Markus Mirz committed
149
150
151
152
153
154
155
156
157
158
159
160
161
162
	}
}

void DP::Ph1::VoltageSource::updateVoltage(Real time) {
	if (mSrcFreq->get() < 0) {
		mIntfVoltage(0,0) = mVoltageRef->get();
	}
	else {
		mIntfVoltage(0,0) = Complex(
			Math::abs(mVoltageRef->get()) * cos(time * 2.*PI*mSrcFreq->get() + Math::phase(mVoltageRef->get())),
			Math::abs(mVoltageRef->get()) * sin(time * 2.*PI*mSrcFreq->get() + Math::phase(mVoltageRef->get())));
	}
}

163
164
165
void DP::Ph1::VoltageSource::mnaPreStep(Real time, Int timeStepCount) {
	updateVoltage(time);
	mnaApplyRightSideVectorStamp(mRightVector);
Markus Mirz's avatar
Markus Mirz committed
166
167
168
169
170
171
172
}

void DP::Ph1::VoltageSource::MnaPreStepHarm::execute(Real time, Int timeStepCount) {
	mVoltageSource.updateVoltage(time);
	mVoltageSource.mnaApplyRightSideVectorStampHarm(mVoltageSource.mRightVector);
}

173
174
void DP::Ph1::VoltageSource::mnaPostStep(Real time, Int timeStepCount, Attribute<Matrix>::Ptr &leftVector) {
	mnaUpdateCurrent(*leftVector);
Markus Mirz's avatar
Markus Mirz committed
175
176
177
178
179
180
181
182
}

void DP::Ph1::VoltageSource::MnaPostStepHarm::execute(Real time, Int timeStepCount) {
	mVoltageSource.mnaUpdateCurrent(*mLeftVectors[0]);
}

void DP::Ph1::VoltageSource::mnaUpdateCurrent(const Matrix& leftVector) {
	for (UInt freq = 0; freq < mNumFreqs; freq++) {
Markus Mirz's avatar
Markus Mirz committed
183
		mIntfCurrent(0,freq) = Math::complexFromVectorElement(leftVector, mVirtualNodes[0]->matrixNodeIndex(), mNumFreqs, freq);
Markus Mirz's avatar
Markus Mirz committed
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
	}
}

void DP::Ph1::VoltageSource::daeResidual(double ttime, const double state[], const double dstate_dt[], double resid[], std::vector<int>& off){
	/* new state vector definintion:
		state[0]=node0_voltage
		state[1]=node1_voltage
		....
		state[n]=noden_voltage
		state[n+1]=component0_voltage
		state[n+2]=component0_inductance (not yet implemented)
		...
		state[m-1]=componentm_voltage
		state[m]=componentm_inductance
	*/

Markus Mirz's avatar
Markus Mirz committed
200
201
    int Pos1 = matrixNodeIndex(0);
    int Pos2 = matrixNodeIndex(1);
Markus Mirz's avatar
Markus Mirz committed
202
203
204
205
206
207
208
209
210
211
212
213
214
215
	int c_offset = off[0]+off[1]; //current offset for component
	int n_offset_1 = c_offset + Pos1 +1;// current offset for first nodal equation
	int n_offset_2 = c_offset + Pos2 +1;// current offset for second nodal equation
	resid[c_offset] = (state[Pos2]-state[Pos1]) - state[c_offset]; // Voltage equation for Resistor
	//resid[++c_offset] = ; //TODO : add inductance equation
	resid[n_offset_1] += mIntfCurrent(0, 0).real();
	resid[n_offset_2] += mIntfCurrent(0, 0).real();
	off[1] += 1;
}

Complex DP::Ph1::VoltageSource::daeInitialize() {
	mIntfVoltage(0,0) = mVoltageRef->get();
	return mVoltageRef->get();
}