Reader.cpp 12.6 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
/** Read CIM files
 *
 * @author Georg Reinke <georg.reinke@rwth-aachen.de>
 * @copyright 2017, Institute for Automation of Complex Power Systems, EONERC
 * @license GNU General Public License (version 3)
 *
 * 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/>.
 *********************************************************************************/

23
24
#include <CIMModel.hpp>
#include <IEC61970.hpp>
Markus Mirz's avatar
Markus Mirz committed
25
26
#include "../Definitions.h"
#include "Reader.h"
27
#include "Components.h"
28
29

using namespace DPsim;
30
using namespace DPsim::CIM;
31
using namespace IEC61970::Base::Core;
32
33
using namespace IEC61970::Base::Domain;
using namespace IEC61970::Base::Equivalents;
34
35
36
37
using namespace IEC61970::Base::Topology;
using namespace IEC61970::Base::Wires;

// TODO is UnitMulitplier actually used/set anywhere?
38
Real Reader::unitValue(Real value, UnitMultiplier mult) {
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
	switch (mult) {
	case UnitMultiplier::p:
		value *= 1e-12;
		break;
	case UnitMultiplier::n:
		value *= 1e-9;
		break;
	case UnitMultiplier::micro:
		value *= 1e-6;
		break;
	case UnitMultiplier::m:
		value *= 1e-3;
		break;
	case UnitMultiplier::c:
		value *= 1e-2;
		break;
	case UnitMultiplier::d:
		value *= 1e-1;
		break;
	case UnitMultiplier::k:
		value *= 1e3;
		break;
	case UnitMultiplier::M:
		value *= 1e6;
		break;
	case UnitMultiplier::G:
		value *= 1e9;
		break;
	case UnitMultiplier::T:
		value *= 1e12;
		break;
	default:
		break;
	}
	return value;
}

76
// TODO: fix error with frequency and angular frequency
77
Reader::Reader(Real systemFrequency, Logger& logger) {
78
	mModel.setDependencyCheckOff();
79
	mNumVoltageSources = 0;
80
	mVoltages = nullptr;
81
	mFrequency = systemFrequency;
82
	mLogger = &logger;
83
84
}

85
Reader::~Reader() {
86
87
	if (mVoltages)
		delete[] mVoltages;
88
89
}

90
ElementPtr Reader::mapACLineSegment(ACLineSegment* line) {
91
	std::vector<Matrix::Index> &nodes = mEqNodeMap.at(line->mRID); // TODO can fail
92
	if (nodes.size() != 2) {
93
		mLogger->Log(LogLevel::WARN) << "ACLineSegment " << line->mRID << " has " << nodes.size() << " terminals, ignoring" << std::endl;
94
95
96
		// TODO better error handling (throw exception?)
		return nullptr;
	}
Steffen Vogel's avatar
Steffen Vogel committed
97

98
99
	Real r = line->r.value;
	Real x = line->x.value;
Steffen Vogel's avatar
Steffen Vogel committed
100

101
	mLogger->Log(LogLevel::INFO) << "Found ACLineSegment " << line->name << " rid=" << line->mRID << " node1=" << nodes[0] << " node2=" << nodes[1]
102
103
104
		<< " R=" << r << " X=" << x << std::endl;

	mLogger->Log(LogLevel::INFO) << "Create RxLine " << line->name << " node1=" << nodes[0] << " node2=" << nodes[1]
105
		<< " R=" << r << " X=" << x << std::endl;
Markus Mirz's avatar
Markus Mirz committed
106
	return std::make_shared<RxLineDP>(line->name, nodes[0], nodes[1], r, x/mFrequency);
107
108
}

109
110
void Reader::mapAsynchronousMachine(AsynchronousMachine* machine) {

111
112
}

113
114
void Reader::mapEnergyConsumer(EnergyConsumer* con) {

115
116
}

117
118
void Reader::mapEquivalentInjection(EquivalentInjection* inj) {

119
120
}

121
ElementPtr Reader::mapExternalNetworkInjection(ExternalNetworkInjection* inj) {
122
	std::vector<Matrix::Index> &nodes = mEqNodeMap.at(inj->mRID);
123
	if (nodes.size() != 1) {
124
		mLogger->Log(LogLevel::ERROR) << "ExternalNetworkInjection " << inj->mRID << " has " << nodes.size() << " terminals, ignoring" << std::endl;
125
126
		return nullptr;
	}
127
	Matrix::Index node = nodes[0];
128
129
	SvVoltage *volt = mVoltages[node-1];
	if (!volt) {
130
		mLogger->Log(LogLevel::ERROR) << "ExternalNetworkInjection " << inj->mRID << " has no associated SvVoltage, ignoring" << std::endl;
131
132
		return nullptr;
	}
133
134
135
	Real voltAbs = volt->v.value;
	Real voltPhase = volt->angle.value;
	Complex initVoltage = std::polar(voltAbs, voltPhase * PI / 180);
136
	mLogger->Log(LogLevel::INFO) << "IdealVoltageSource " << inj->name << " rid=" << inj->mRID << " node1=" << node 
137
138
		<< " V=" << voltAbs << "<" << voltPhase << std::endl;
	
Markus Mirz's avatar
Markus Mirz committed
139
	return std::make_shared<IdealVoltageSource>(inj->name, node, 0, initVoltage);
140
141
}

142
ElementPtr Reader::mapPowerTransformer(PowerTransformer* trans) {
143
	std::vector<Matrix::Index> &nodes = mEqNodeMap.at(trans->mRID);
144
	if (nodes.size() != trans->PowerTransformerEnd.size()) {
145
		mLogger->Log(LogLevel::WARN) << "PowerTransformer " << trans->mRID << " has differing number of terminals and windings, ignoring" << std::endl;
146
147
148
149
		return nullptr;
	}
	if (nodes.size() != 2) {
		// TODO three windings also possible
150
		mLogger->Log(LogLevel::WARN) << "PowerTransformer " << trans->mRID << " has " << nodes.size() << "terminals; ignoring" << std::endl;
151
152
		return nullptr;
	}
153

154
	mLogger->Log(LogLevel::INFO) << "Found PowerTransformer " << trans->name << " rid=" << trans->mRID
155
		<< " node1=" << nodes[0] << " node2=" << nodes[1] << std::endl;
156

157
158
	Real voltageNode1 = 0;
	Real voltageNode2 = 0;
159
160
	for (PowerTransformerEnd *end : trans->PowerTransformerEnd) {
		if (end->endNumber == 1) {
161
			mLogger->Log(LogLevel::INFO) << "    PowerTransformerEnd_1 " << end->name
162
163
				<< " Vrated=" << end->ratedU.value << " R=" << end->r.value << " X=" << end->x.value << std::endl;
			voltageNode1 = end->ratedU.value;
164
		}
165
166
167
168
169
		if (end->endNumber == 2) {
			mLogger->Log(LogLevel::INFO) << "    PowerTransformerEnd_2 " << end->name
				<< " Vrated=" << end->ratedU.value << " R=" << end->r.value << " X=" << end->x.value << std::endl;
			voltageNode2 = end->ratedU.value;
		}
170
	}
171

172
	if (voltageNode1 != 0 && voltageNode2 != 0) {
Markus Mirz's avatar
Markus Mirz committed
173
174
		Real ratioAbs = voltageNode1 / voltageNode2;
		mLogger->Log(LogLevel::INFO) << "    Calculated ratio=" << ratioAbs << std::endl;
175
176
		mLogger->Log(LogLevel::INFO) << "Create PowerTransformer " << trans->name
			<< " node1=" << nodes[0] << " node2=" << nodes[1]
Markus Mirz's avatar
Markus Mirz committed
177
178
			<< " ratio=" << ratioAbs << "<0" << std::endl;
		Real ratioPhase = 0;
179
		return std::make_shared<IdealTransformerDP>(trans->name, nodes[0], nodes[1], ratioAbs, ratioPhase);
180

181
	}
182
	mLogger->Log(LogLevel::WARN) << "PowerTransformer " << trans->mRID << " has no primary End; ignoring" << std::endl;
183
184
185
	return nullptr;
}

Markus Mirz's avatar
Markus Mirz committed
186
// TODO: don't use SvVoltage, but map to a SynchronGeneratorDP instead
187
ElementPtr Reader::mapSynchronousMachine(SynchronousMachine* machine) {
188
	std::vector<Matrix::Index> &nodes = mEqNodeMap.at(machine->mRID);
189
	if (nodes.size() != 1) {
190
191
		// TODO: check with the model if this assumption (only 1 terminal) is always true
		mLogger->Log(LogLevel::WARN) << "SynchronousMachine " << machine->mRID << " has " << nodes.size() << " terminals, ignoring" << std::endl;
192
193
		return nullptr;
	}
194
	Matrix::Index node = nodes[0];
195
196
	SvVoltage *volt = mVoltages[node-1];
	if (!volt) {
197
		mLogger->Log(LogLevel::WARN) << "SynchronousMachine " << machine->mRID << " has no associated SvVoltage, ignoring" << std::endl;
198
199
		return nullptr;
	}
200

201
202
203
204
	Real voltAbs = volt->v.value;
	Real voltPhase = volt->angle.value;
	Complex initVoltage = std::polar(voltAbs, voltPhase * PI / 180);
	
205
	// Apply unit multipliers according to CGMES convetions.
206
	volt->v.value = Reader::unitValue(volt->v.value, UnitMultiplier::k);
207

208
	// TODO is it appropiate to use this resistance here
209
	mLogger->Log(LogLevel::INFO) << "Create IdealVoltageSource " << machine->name << " node=" << node
210
		<< " V=" << voltAbs << "<" << voltPhase << std::endl;
Markus Mirz's avatar
Markus Mirz committed
211
	return std::make_shared<IdealVoltageSource>(machine->name, node, 0, initVoltage);
212
213
}

214
ElementPtr Reader::mapComponent(BaseClass* obj) {
215
	if (ACLineSegment *line = dynamic_cast<ACLineSegment*>(obj))
216
		return mapACLineSegment(line);
217
218
	if (ExternalNetworkInjection *inj = dynamic_cast<ExternalNetworkInjection*>(obj))
		return mapExternalNetworkInjection(inj);
219
220
	if (PowerTransformer *trans = dynamic_cast<PowerTransformer*>(obj))
		return mapPowerTransformer(trans);
221
222
	if (SynchronousMachine *syncMachine = dynamic_cast<SynchronousMachine*>(obj))
		return mapSynchronousMachine(syncMachine);
223
224
225
	return nullptr;
}

226
ElementPtr Reader::newPQLoad(String rid, String name) {
227
	std::vector<Matrix::Index> &nodes = mEqNodeMap.at(rid);
228
	if (nodes.size() != 1) {
229
		mLogger->Log(LogLevel::WARN) << rid << " has " << nodes.size() << " terminals; ignoring" << std::endl;
230
231
232
233
		return nullptr;
	}
	auto search = mPowerFlows.find(rid);
	if (search == mPowerFlows.end()) {
234
		mLogger->Log(LogLevel::WARN) << rid << " has no associated SvPowerFlow, ignoring" << std::endl;
235
236
237
		return nullptr;
	}
	SvPowerFlow* flow = search->second;
238
	Matrix::Index node = nodes[0];
239
240
	SvVoltage *volt = mVoltages[node-1];
	if (!volt) {
241
		mLogger->Log(LogLevel::WARN) << rid << " has no associated SvVoltage, ignoring" << std::endl;
242
243
		return nullptr;
	}
244

245
246
	mLogger->Log(LogLevel::INFO) << "Found EnergyConsumer " << name << " rid=" << rid << " node="
		<< node << " P=" << flow->p.value << " Q=" << flow->q.value
247
248
249
		<< " V=" << volt->v.value << "<" << volt->angle.value << std::endl;

	// Apply unit multipliers according to CGMES convetions.
250
251
252
253
	flow->p.value = Reader::unitValue(flow->p.value, UnitMultiplier::M);
	flow->q.value = Reader::unitValue(flow->q.value, UnitMultiplier::M);
	volt->v.value = Reader::unitValue(volt->v.value, UnitMultiplier::k);

254
255
256
	mLogger->Log(LogLevel::INFO) << "Create PQLoad " << name << " node="
		<< node << " P=" << flow->p.value << " Q=" << flow->q.value
		<< " V=" << volt->v.value << "<" << volt->angle.value << std::endl;
257
	return std::make_shared<PQLoadDP>(name, node, flow->p.value, flow->q.value, volt->v.value, volt->angle.value*PI/180);
258
259
}

260
bool Reader::addFile(String filename) {
261
262
263
	return mModel.addCIMFile(filename);
}

264
void Reader::parseFiles() {
265
266
267
268
	mModel.parseFiles();
	// First, go through all topological nodes and collect them in a list.
	// Since all nodes have references to the equipment connected to them (via Terminals), but not
	// the other way around (which we need for instantiating the components), we collect that information here as well.
269
	mLogger->Log(LogLevel::INFO) << "#### List of topological nodes and associated terminals ####" << std::endl;
270
271
272
	for (BaseClass* obj : mModel.Objects) {
		TopologicalNode* topNode = dynamic_cast<TopologicalNode*>(obj);
		if (topNode) {
273
			mLogger->Log(LogLevel::INFO) << "TopologicalNode " << mTopNodes.size()+1 << " rid=" << topNode->mRID << " Terminals:" << std::endl;
274
			mTopNodes[topNode->mRID] = (Matrix::Index) mTopNodes.size()+1;
275
			for (Terminal* term : topNode->Terminal) {
276
				mLogger->Log(LogLevel::INFO) << "    " << term->mRID << std::endl;
277
278
				ConductingEquipment *eq = term->ConductingEquipment;
				if (!eq) {
279
					mLogger->Log(LogLevel::WARN) << "Terminal " << term->mRID << " has no Conducting Equipment, ignoring!" << std::endl;
280
				} else {
281
					mLogger->Log(LogLevel::INFO) << "    eq " << eq->mRID << " sequenceNumber " << term->sequenceNumber << std::endl;
282
283
					std::vector<Matrix::Index> &nodesVec = mEqNodeMap[eq->mRID];
					if (nodesVec.size() < (unsigned) term->sequenceNumber) {
284
						nodesVec.resize(term->sequenceNumber);
285
					}
286
					nodesVec[term->sequenceNumber-1] = (Matrix::Index) mTopNodes.size();
287
288
289
290
				}
			}
		}
	}
291
292
293
	// Collect voltage state variables associated to nodes that are used
	// for various components.
	mVoltages = new SvVoltage*[mTopNodes.size()];
294
	mLogger->Log(LogLevel::INFO) << "#### List of node voltages from power flow calculation ####" << std::endl;
295
	for (BaseClass* obj : mModel.Objects) {
296
		if (SvVoltage* volt = dynamic_cast<SvVoltage*>(obj)) {
297
298
			TopologicalNode* node = volt->TopologicalNode;
			if (!node) {
299
				mLogger->Log(LogLevel::WARN) << "SvVoltage references missing Topological Node, ignoring" << std::endl;
300
301
302
303
				continue;
			}
			auto search = mTopNodes.find(node->mRID);
			if (search == mTopNodes.end()) {
304
				mLogger->Log(LogLevel::WARN) << "SvVoltage references Topological Node " << node->mRID << " missing from mTopNodes, ignoring" << std::endl;
305
306
307
				continue;
			}
			mVoltages[search->second-1] = volt;
308
			mLogger->Log(LogLevel::INFO) << "Node " << search->second << ": " << volt->v.value << "<" << volt->angle.value << std::endl;
309
310
311
312
		} else if (SvPowerFlow* flow = dynamic_cast<SvPowerFlow*>(obj)) {
			// TODO could there be more than one power flow per equipment?
			Terminal* term = flow->Terminal;
			mPowerFlows[term->ConductingEquipment->mRID] = flow;
313
314
		}
	}
315
	mLogger->Log(LogLevel::INFO) << "#### Create new components ####" << std::endl;
316
	for (BaseClass* obj : mModel.Objects) {
Markus Mirz's avatar
Markus Mirz committed
317
		ElementPtr comp = mapComponent(obj);
318
		if (comp)
319
			mComponents.push_back(comp);
320
	}
321
322
}

323
ElementList& Reader::getComponents() {
324
325
326
	return mComponents;
}

327
Matrix::Index Reader::mapTopologicalNode(String mrid) {
328
329
330
	auto search = mTopNodes.find(mrid);
	if (search == mTopNodes.end())
		return -1;
Steffen Vogel's avatar
Steffen Vogel committed
331

332
333
334
	return search->second;
}

335
int Reader::getNumVoltageSources() {
336
	return mNumVoltageSources;
337
}