Power System Platform  2026w10a-beta
Loading...
Searching...
No Matches
PowerFlow.cpp
1/*
2 * Copyright (C) 2017 Thales Lima Oliveira <thales@ufu.br>
3 *
4 * This program is free software; you can redistribute it and/or modify
5 * it under the terms of the GNU General Public License as published by
6 * the Free Software Foundation; either version 2 of the License, or
7 * any later version.
8 *
9 * This program is distributed in the hope that it will be useful,
10 * but WITHOUT ANY WARRANTY; without even the implied warranty of
11 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12 * GNU General Public License for more details.
13 *
14 * You should have received a copy of the GNU General Public License
15 * along with this program. If not, see <https://www.gnu.org/licenses/>.
16 */
17
18#include "PowerFlow.h"
19
20PowerFlow::PowerFlow() : ElectricCalculation() {}
21PowerFlow::PowerFlow(std::vector<Element*> elementList) : ElectricCalculation() { GetElementsFromList(elementList); }
22PowerFlow::~PowerFlow() {}
23bool PowerFlow::InitPowerFlow(std::vector<BusType>& busType,
24 std::vector<std::complex<double> >& voltage,
25 std::vector<std::complex<double> >& power,
26 std::vector<std::complex<double> >& loadPower,
27 std::vector<ReactiveLimits>& reactiveLimit,
28 double systemPowerBase,
29 double initAngle)
30{
31 double radInitAngle = wxDegToRad(initAngle);
32
33 // Calculate EMT Elements admittance to add to the Ybus.
34 if (!CalculateEMTElementsAdmittance(systemPowerBase, m_errorMsg)) return false;
35
36 // Calculate the Ybus.
37 if (!GetYBus(m_yBus, systemPowerBase)) {
38 m_errorMsg = _("No buses found on the system.");
39 return false;
40 }
41
42 // Calculate EMT elements power for 1.0 p.u. voltage
43 if (!CalculateEMTElementsPower(systemPowerBase, m_errorMsg, false)) return false;
44
45 // Number of buses in the system.
46 //m_numberOfBuses = static_cast<int>(m_busList.size());
47 m_numberOfBuses = static_cast<int>(m_yBus.size());
48
49 busType.clear();
50 voltage.clear();
51 power.clear();
52 loadPower.clear();
53 reactiveLimit.clear();
54
55 reactiveLimit.resize(m_numberOfBuses);
56
57 int busNumber = 0;
58 Bus* slackBus = nullptr;
59 for (auto itb = m_busList.begin(); itb != m_busList.end(); itb++) {
60 Bus* bus = *itb;
61 BusElectricalData data = bus->GetElectricalData();
62 if (data.isConnected) {
63 // Fill the bus type
64 if (data.slackBus) {
65 busType.push_back(BUS_SLACK);
66 slackBus = bus;
67 }
68 // If the bus have controlled voltage, check if at least one synchronous machine is connected, then set the
69 // bus type.
70 else if (data.isVoltageControlled) {
71 bool hasSyncMachine = false;
72 // Synchronous generator
73 for (auto itsg = m_syncGeneratorList.begin(); itsg != m_syncGeneratorList.end(); itsg++) {
74 SyncGenerator* syncGenerator = *itsg;
75 if (bus == syncGenerator->GetParentList()[0] && syncGenerator->IsOnline()) hasSyncMachine = true;
76 }
77 // Synchronous motor
78 for (auto itsm = m_syncMotorList.begin(); itsm != m_syncMotorList.end(); itsm++) {
79 SyncMotor* syncMotor = *itsm;
80 if (bus == syncMotor->GetParentList()[0] && syncMotor->IsOnline()) hasSyncMachine = true;
81 }
82 if (hasSyncMachine)
83 busType.push_back(BUS_PV);
84 else
85 busType.push_back(BUS_PQ);
86 }
87 else
88 busType.push_back(BUS_PQ);
89
90 // Fill the voltages array
91 double v = 1.0;
92 double t = 0.0;
93 if (data.slackBus) {
94 v = data.controlledVoltage;
95 t = radInitAngle;
96 }
97 else if (busType[busNumber] == BUS_PV) {
98 v = data.controlledVoltage;
99 t = std::arg(data.voltage);
100 }
101 else {
102 v = std::abs(data.voltage);
103 if (v <= 0.1) v = 1.0; // Avoid very low voltage magnitude that can cause convergence problems.
104 t = std::arg(data.voltage);
105 }
106 voltage.push_back(std::complex<double>(v * std::cos(t), v * std::sin(t)));
107
108 //if (data.isVoltageControlled && busType[busNumber] != BUS_PQ) {
109 // voltage.push_back(std::complex<double>(data.controlledVoltage * std::cos(radInitAngle),
110 // data.controlledVoltage * std::sin(radInitAngle)));
111 //}
112 //else {
113 // voltage.push_back(std::complex<double>(std::cos(radInitAngle), std::sin(radInitAngle)));
114 //}
115
116 // Fill the power array
117 power.push_back(std::complex<double>(0.0, 0.0)); // Initial value
118 loadPower.push_back(std::complex<double>(0.0, 0.0));
119
120 // Synchronous generator
121 for (auto itsg = m_syncGeneratorList.begin(); itsg != m_syncGeneratorList.end(); itsg++) {
122 SyncGenerator* syncGenerator = *itsg;
123 if (syncGenerator->IsOnline()) {
124 if (bus == syncGenerator->GetParentList()[0]) {
125 SyncGeneratorElectricalData childData = syncGenerator->GetPUElectricalData(systemPowerBase);
126 power[busNumber] += std::complex<double>(childData.activePower, childData.reactivePower);
127
128 if (busType[busNumber] == BUS_PV) {
129 if (childData.haveMaxReactive && reactiveLimit[busNumber].maxLimitType != RL_UNLIMITED_SOURCE) {
130 reactiveLimit[busNumber].maxLimitType = RL_LIMITED;
131 reactiveLimit[busNumber].maxLimit += childData.maxReactive;
132 }
133 else if (!childData.haveMaxReactive)
134 reactiveLimit[busNumber].maxLimitType = RL_UNLIMITED_SOURCE;
135
136 if (childData.haveMinReactive && reactiveLimit[busNumber].minLimitType != RL_UNLIMITED_SOURCE) {
137 reactiveLimit[busNumber].minLimitType = RL_LIMITED;
138 reactiveLimit[busNumber].minLimit += childData.minReactive;
139 }
140 else if (!childData.haveMinReactive)
141 reactiveLimit[busNumber].minLimitType = RL_UNLIMITED_SOURCE;
142 }
143 }
144 }
145 }
146 // Synchronous motor
147 for (auto itsm = m_syncMotorList.begin(); itsm != m_syncMotorList.end(); itsm++) {
148 SyncMotor* syncMotor = *itsm;
149 if (syncMotor->IsOnline()) {
150 if (bus == syncMotor->GetParentList()[0]) {
151 SyncMotorElectricalData childData = syncMotor->GetPUElectricalData(systemPowerBase);
152 power[busNumber] += std::complex<double>(-childData.activePower, childData.reactivePower);
153 loadPower[busNumber] += std::complex<double>(-childData.activePower, 0.0);
154
155 if (busType[busNumber] == BUS_PV) {
156 if (childData.haveMaxReactive && reactiveLimit[busNumber].maxLimitType != RL_UNLIMITED_SOURCE) {
157 reactiveLimit[busNumber].maxLimitType = RL_LIMITED;
158 reactiveLimit[busNumber].maxLimit += childData.maxReactive;
159 }
160 else if (!childData.haveMaxReactive)
161 reactiveLimit[busNumber].maxLimitType = RL_UNLIMITED_SOURCE;
162
163 if (childData.haveMinReactive && reactiveLimit[busNumber].minLimitType != RL_UNLIMITED_SOURCE) {
164 reactiveLimit[busNumber].minLimitType = RL_LIMITED;
165 reactiveLimit[busNumber].minLimit += childData.minReactive;
166 }
167 else if (!childData.haveMinReactive)
168 reactiveLimit[busNumber].minLimitType = RL_UNLIMITED_SOURCE;
169 }
170 }
171 }
172 }
173 // Load
174 for (auto itl = m_loadList.begin(); itl != m_loadList.end(); itl++) {
175 Load* load = *itl;
176 if (load->IsOnline()) {
177 if (bus == load->GetParentList()[0]) {
178 LoadElectricalData childData = load->GetPUElectricalData(systemPowerBase);
179 if (childData.loadType == CONST_POWER) {
180 power[busNumber] += std::complex<double>(-childData.activePower, -childData.reactivePower);
181 loadPower[busNumber] += std::complex<double>(-childData.activePower, -childData.reactivePower);
182 }
183 }
184 }
185 }
186
187 // Induction motor
188 for (auto itim = m_indMotorList.begin(); itim != m_indMotorList.end(); itim++) {
189 IndMotor* indMotor = *itim;
190 if (indMotor->IsOnline()) {
191 if (bus == indMotor->GetParentList()[0]) {
192 IndMotorElectricalData childData = indMotor->GetPUElectricalData(systemPowerBase);
193 double reactivePower = childData.reactivePower;
194
195 if (childData.calcQInPowerFlow) {
196 indMotor->InitPowerFlowMotor(systemPowerBase, data.number);
197 if (!indMotor->CalculateReactivePower(std::abs(voltage[childData.busNum]))) {
198 m_errorMsg = _("It was not possible to solve the induction motors.");
199 return false;
200 }
201 reactivePower = indMotor->GetElectricalData().qValue;
202 }
203
204 power[busNumber] += std::complex<double>(-childData.activePower, -reactivePower);
205 loadPower[busNumber] += std::complex<double>(-childData.activePower, -reactivePower);
206 }
207 }
208 }
209
210 // EMTElements
211 for (EMTElement* emtElement : m_emtElementList) {
212 if (emtElement->IsOnline() && !emtElement->GetParentList().empty()) {
213 if (bus == emtElement->GetParentList()[0]) {
214 EMTElementData childData = emtElement->GetEMTElementData();
215 power[busNumber] -= childData.power;
216 loadPower[busNumber] -= childData.power;
217 }
218 }
219 }
220
221 busNumber++;
222 }
223 }
224
225 // Check if have slack bus and if have generation on the slack bus
226 bool haveSlackBus = false;
227 bool slackBusHaveGeneration = false;
228 for (unsigned int i = 0; i < busType.size(); i++) {
229 if (busType[i] == BUS_SLACK) {
230 for (auto itsg = m_syncGeneratorList.begin(); itsg != m_syncGeneratorList.end(); itsg++) {
231 SyncGenerator* syncGenerator = *itsg;
232 if (syncGenerator->IsOnline() && slackBus == syncGenerator->GetParentList()[0]) slackBusHaveGeneration = true;
233 }
234 haveSlackBus = true;
235 }
236 }
237 if (!haveSlackBus) {
238 m_errorMsg = _("There is no slack bus on the system.");
239 return false;
240 }
241 if (!slackBusHaveGeneration) {
242 m_errorMsg = _("The slack bus don't have generation.");
243 return false;
244 }
245
246 return true;
247}
248
249bool PowerFlow::RunGaussSeidel(double systemPowerBase,
250 int maxIteration,
251 double error,
252 double initAngle,
253 double accFactor)
254{
255 std::vector<BusType> busType; // Bus type
256 std::vector<std::complex<double> > voltage; // Voltage of buses
257 std::vector<std::complex<double> > power; // Injected power
258 std::vector<std::complex<double> > loadPower; // Only the load power
259 std::vector<ReactiveLimits> reactiveLimit; // Limit of reactive power on PV buses
260
261 if (!InitPowerFlow(busType, voltage, power, loadPower, reactiveLimit, systemPowerBase, initAngle)) return false;
262
263 // Gauss-Seidel method
264 std::vector<std::complex<double> > oldVoltage; // Old voltage array.
265 oldVoltage.resize(voltage.size());
266
267 auto oldBusType = busType;
268
269 int iteration = 0; // Current itaration number.
270 double emtPowerError = 1e3; // EMT Elements power error.
271
272 while (true) {
273 // Reach the max number of iterations.
274 if (iteration >= maxIteration) {
275 m_errorMsg = _("The maximum number of iterations was reached.");
276 return false;
277 }
278
279 // Calculate induction motor reactive power
280 if (!CalculateMotorsReactivePower(voltage, power)) {
281 m_errorMsg = _("It was not possible to solve the induction motors.");
282 return false;
283 }
284
285 // Update the old voltage array to current iteration values.
286 for (int i = 0; i < m_numberOfBuses; i++) oldVoltage[i] = voltage[i];
287
288 double iterationError = GaussSeidel(busType, voltage, oldVoltage, power, accFactor);
289
290 if (HasInvalidValue(voltage)) {
291 m_errorMsg = _("The power flow solution has invalid voltage values.");
292 ResetVoltages();
293 return false;
294 }
295 if (HasInvalidValue(power)) {
296 m_errorMsg = _("The power flow solution has invalid power values.");
297 ResetVoltages();
298 return false;
299 }
300
301 if (iterationError < error) {
302 // Calculate power error in EMT elements.
303 if (!m_emtElementList.empty()) {
304 double oldError = emtPowerError;
305 emtPowerError = CalculateEMTPowerError(voltage, power, systemPowerBase, m_errorMsg);
306 //wxMessageBox(wxString::Format(wxT("It %d, EMT Power Error: %e (base = %e)"), iteration, emtPowerError, error));
307 if (abs(emtPowerError - oldError) < 1e-12) { // Without change in error the convergence is impossible
308 m_errorMsg = _("Impossible to reach a convergence with the current Electromagnetic Transient Elements.");
309 return false;
310 }
311 }
312 else emtPowerError = 0.0;
313
314 if (!CheckReactiveLimits(busType, reactiveLimit, power, loadPower) && emtPowerError < error) break;
315 }
316
317 iteration++;
318 }
319 m_iterations = iteration;
320
321 // Adjust the power array.
322 for (int i = 0; i < m_numberOfBuses; i++) {
323 std::complex<double> sBus = std::complex<double>(0.0, 0.0);
324 for (int j = 0; j < m_numberOfBuses; j++) sBus += voltage[i] * std::conj(voltage[j]) * std::conj(m_yBus[i][j]);
325 power[i] = sBus;
326 }
327
328
329 UpdateElementsPowerFlow(voltage, power, oldBusType, reactiveLimit, systemPowerBase);
330
331 return true;
332}
333
334bool PowerFlow::HasInvalidValue(const std::vector< std::complex<double> >& value)
335{
336 for (size_t i = 0; i < value.size(); i++)
337 {
338 if (!std::isfinite(value[i].real()) ||
339 !std::isfinite(value[i].imag()))
340 {
341 return true;
342 }
343 }
344 return false;
345}
346
347bool PowerFlow::RunNewtonRaphson(double systemPowerBase,
348 int maxIteration,
349 double error,
350 double initAngle,
351 double inertia)
352{
353 std::vector<BusType> busType; // Bus type
354 std::vector<std::complex<double> > voltage; // Voltage of buses
355 std::vector<std::complex<double> > power; // Injected power
356 std::vector<std::complex<double> > loadPower; // Only the load power
357 std::vector<ReactiveLimits> reactiveLimit; // Limit of reactive power on PV buses
358
359 if (!InitPowerFlow(busType, voltage, power, loadPower, reactiveLimit, systemPowerBase, initAngle)) return false;
360 auto oldBusType = busType;
361
362 // Newton-Raphson method
363 int numPQ = 0; // Number of PQ buses
364 int numPV = 0; // Number of PV buses
365 GetNumPVPQ(busType, numPQ, numPV);
366
367 // DeltaP and DeltaQ array
368 std::vector<double> dPdQ;
369 dPdQ.resize(numPV + 2 * numPQ, 0.0);
370
371 int iteration = 0; // Current iteration number.
372 while (true) {
373 // Reach the max number of iterations.
374 if (iteration >= maxIteration) {
375 m_errorMsg = _("The maximum number of iterations was reached.");
376 return false;
377 }
378
379 // Calculate induction motor reactive power
380 if (!CalculateMotorsReactivePower(voltage, power)) {
381 m_errorMsg = _("It was not possible to solve the induction motors.");
382 return false;
383 }
384
385 // Calculate dPdQ array
386
387 // Fill it with zeros
388 std::fill(dPdQ.begin(), dPdQ.end(), 0.0);
389
390 int indexDP = 0;
391 int indexDQ = numPQ + numPV;
392 for (int i = 0; i < m_numberOfBuses; i++) {
393 if (busType[i] != BUS_SLACK) {
394 for (int j = 0; j < m_numberOfBuses; j++) {
395 // PV ou PQ bus
396 std::complex<double> sInj = std::conj(m_yBus[i][j]) * voltage[i] * std::conj(voltage[j]);
397 dPdQ[indexDP] += sInj.real();
398
399 // PQ bus
400 if (busType[i] == BUS_PQ) dPdQ[indexDQ] += sInj.imag();
401 }
402
403 // PQ or PV bus
404 dPdQ[indexDP] = power[i].real() - dPdQ[indexDP];
405 indexDP++;
406
407 // PQ bus
408 if (busType[i] == BUS_PQ) {
409 dPdQ[indexDQ] = power[i].imag() - dPdQ[indexDQ];
410 indexDQ++;
411 }
412 }
413 }
414
415 // Calculate the iteration error
416 double iterationError = 0.0;
417 for (unsigned int i = 0; i < dPdQ.size(); ++i) {
418 if (iterationError < std::abs(dPdQ[i])) iterationError = std::abs(dPdQ[i]);
419 }
420
421 // Check if the iteration error is less than tolerance, also check if any reactive limit was reached.
422 // If any reactive limit was reached, change the bus type.
423 if (iterationError < error) {
424 for (int i = 0; i < m_numberOfBuses; i++) {
425 std::complex<double> sBus(0.0, 0.0);
426 for (int j = 0; j < m_numberOfBuses; j++)
427 sBus += voltage[i] * std::conj(voltage[j]) * std::conj(m_yBus[i][j]);
428 power[i] = sBus;
429 }
430
431 double emtPowerError = CalculateEMTPowerError(voltage, power, systemPowerBase, m_errorMsg);
432
433 if (!CheckReactiveLimits(busType, reactiveLimit, power, loadPower) && emtPowerError < error)
434 break;
435 else {
436 GetNumPVPQ(busType, numPQ, numPV);
437 dPdQ.assign(numPV + 2 * numPQ, 0.0);
438 }
439 }
440
441 NewtonRaphson(busType, voltage, power, numPV, numPQ, dPdQ, inertia);
442
443 if (HasInvalidValue(voltage)) {
444 m_errorMsg = _("The power flow solution has invalid voltage values.");
445 ResetVoltages();
446 return false;
447 }
448 if (HasInvalidValue(power)) {
449 m_errorMsg = _("The power flow solution has invalid power values.");
450 ResetVoltages();
451 return false;
452 }
453
454 iteration++;
455 }
456 m_iterations = iteration;
457
458
459
460 // Adjust the power array.
461 for (int i = 0; i < m_numberOfBuses; i++) {
462 std::complex<double> sBus = std::complex<double>(0.0, 0.0);
463 for (int j = 0; j < m_numberOfBuses; j++) sBus += voltage[i] * std::conj(voltage[j]) * std::conj(m_yBus[i][j]);
464 power[i] = sBus;
465 }
466
467 UpdateElementsPowerFlow(voltage, power, oldBusType, reactiveLimit, systemPowerBase);
468
469 return true;
470}
471
472bool PowerFlow::RunGaussNewton(double systemPowerBase,
473 int maxIteration,
474 double error,
475 double initAngle,
476 double accFactor,
477 double gaussTol,
478 double inertia)
479{
480 std::vector<BusType> busType; // Bus type
481 std::vector<std::complex<double> > voltage; // Voltage of buses
482 std::vector<std::complex<double> > power; // Injected power
483 std::vector<std::complex<double> > loadPower; // Only the load power
484 std::vector<ReactiveLimits> reactiveLimit; // Limit of reactive power on PV buses
485
486 if (!InitPowerFlow(busType, voltage, power, loadPower, reactiveLimit, systemPowerBase, initAngle)) return false;
487
488 // Gauss-Seidel method
489 std::vector<std::complex<double> > oldVoltage; // Old voltage array.
490 oldVoltage.resize(voltage.size());
491
492 auto oldBusType = busType;
493
494 int iteration = 0; // Current itaration number.
495
496 while (true) {
497 // Reach the max number of iterations.
498 if (iteration >= maxIteration) {
499 m_errorMsg = _("The maximum number of iterations was reached.");
500 return false;
501 }
502
503 // Calculate induction motor reactive power
504 if (!CalculateMotorsReactivePower(voltage, power)) {
505 m_errorMsg = _("It was not possible to solve the induction motors.");
506 return false;
507 }
508
509 // Update the old voltage array to current iteration values.
510 for (int i = 0; i < m_numberOfBuses; i++) oldVoltage[i] = voltage[i];
511
512 double iterationError = GaussSeidel(busType, voltage, oldVoltage, power, accFactor);
513
514 if (iterationError < gaussTol) break;
515
516 iteration++;
517 }
518
519 // Newton-Raphson method
520 int numPQ = 0; // Number of PQ buses
521 int numPV = 0; // Number of PV buses
522 GetNumPVPQ(busType, numPQ, numPV);
523
524 // DeltaP and DeltaQ array
525 std::vector<double> dPdQ;
526 dPdQ.resize(numPV + 2 * numPQ, 0.0);
527
528 while (true) {
529 // Reach the max number of iterations.
530 if (iteration >= maxIteration) {
531 m_errorMsg = _("The maximum number of iterations was reached.");
532 return false;
533 }
534
535 // Calculate induction motor reactive power
536 if (!CalculateMotorsReactivePower(voltage, power)) {
537 m_errorMsg = _("It was not possible to solve the induction motors.");
538 return false;
539 }
540
541 // Calculate dPdQ array
542
543 // Fill it with zeros
544 std::fill(dPdQ.begin(), dPdQ.end(), 0.0);
545
546 int indexDP = 0;
547 int indexDQ = numPQ + numPV;
548 for (int i = 0; i < m_numberOfBuses; i++) {
549 if (busType[i] != BUS_SLACK) {
550 for (int j = 0; j < m_numberOfBuses; j++) {
551 // PV ou PQ bus
552 std::complex<double> sInj = std::conj(m_yBus[i][j]) * voltage[i] * std::conj(voltage[j]);
553 dPdQ[indexDP] += sInj.real();
554
555 // PQ bus
556 if (busType[i] == BUS_PQ) dPdQ[indexDQ] += sInj.imag();
557 }
558
559 // PQ or PV bus
560 dPdQ[indexDP] = power[i].real() - dPdQ[indexDP];
561 indexDP++;
562
563 // PQ bus
564 if (busType[i] == BUS_PQ) {
565 dPdQ[indexDQ] = power[i].imag() - dPdQ[indexDQ];
566 indexDQ++;
567 }
568 }
569 }
570
571 // Calculate the iteration error
572 double iterationError = 0.0;
573 for (unsigned int i = 0; i < dPdQ.size(); ++i) {
574 if (iterationError < std::abs(dPdQ[i])) iterationError = std::abs(dPdQ[i]);
575 }
576
577 // Check if the iteration error is less than tolerance, also check if any reactive limit was reached.
578 // If any reactive limit was reached, change the bus type.
579 if (iterationError < error) {
580 double emtPowerError = CalculateEMTPowerError(voltage, power, systemPowerBase, m_errorMsg);
581
582 if (!CheckReactiveLimits(busType, reactiveLimit, power, loadPower) && emtPowerError < error)
583 break;
584 else {
585 GetNumPVPQ(busType, numPQ, numPV);
586 dPdQ.clear();
587 dPdQ.resize(numPV + 2 * numPQ, 0.0);
588 }
589 }
590
591 NewtonRaphson(busType, voltage, power, numPV, numPQ, dPdQ, inertia);
592
593 iteration++;
594 }
595 m_iterations = iteration;
596
597 // Adjust the power array.
598 for (int i = 0; i < m_numberOfBuses; i++) {
599 std::complex<double> sBus = std::complex<double>(0.0, 0.0);
600 for (int j = 0; j < m_numberOfBuses; j++) sBus += voltage[i] * std::conj(voltage[j]) * std::conj(m_yBus[i][j]);
601 power[i] = sBus;
602 }
603
604 UpdateElementsPowerFlow(voltage, power, oldBusType, reactiveLimit, systemPowerBase);
605
606 return true;
607}
608
609void PowerFlow::ResetVoltages()
610{
611 for (auto* bus : m_busList) {
612 BusElectricalData data = bus->GetElectricalData();
613 if (bus->IsOnline()) {
614 data.voltage = std::complex<double>(1.0, 0.0);
615 if (data.isVoltageControlled) {
616 data.voltage = std::complex<double>(data.controlledVoltage, 0.0);
617 }
618 }
619 else
620 data.voltage = std::complex<double>(0.0, 0.0);
621 data.harmonicOrder.clear();
622 data.harmonicVoltage.clear();
623 data.thd = 0.0;
624 bus->SetElectricalData(data);
625 }
626}
627
628void PowerFlow::GetNumPVPQ(std::vector<BusType> busType, int& numPQ, int& numPV)
629{
630 numPQ = 0;
631 numPV = 0;
632 for (auto it = busType.begin(), itEnd = busType.end(); it != itEnd; ++it) {
633 if (*it == BUS_PQ)
634 numPQ++;
635 else if (*it == BUS_PV)
636 numPV++;
637 }
638}
639
640double PowerFlow::GaussSeidel(std::vector<BusType> busType,
641 std::vector<std::complex<double> >& voltage,
642 std::vector<std::complex<double> > oldVoltage,
643 std::vector<std::complex<double> >& power,
644 double accFactor)
645{
646 double error = 0.0;
647
648 for (int i = 0; i < m_numberOfBuses; i++) {
649 if (busType[i] == BUS_PQ) {
650 std::complex<double> yeSum(0.0, 0.0);
651 for (int k = 0; k < m_numberOfBuses; k++) {
652 if (i != k) {
653 // Sum { Y[i,k] * E[k] } | k = 1->n; k diff i
654 yeSum += m_yBus[i][k] * voltage[k];
655 }
656 }
657
658 // E[i] = (1/Y[i,i])*((P[i]-jQ[i])/E*[i] - Sum { Y[i,k] * E[k] (k diff i) })
659 std::complex<double> newVolt = (1.0 / m_yBus[i][i]) * (std::conj(power[i]) / std::conj(voltage[i]) - yeSum);
660
661 // Apply the acceleration factor.
662 newVolt = std::complex<double>(accFactor * (newVolt.real() - voltage[i].real()) + voltage[i].real(),
663 accFactor * (newVolt.imag() - voltage[i].imag()) + voltage[i].imag());
664
665 voltage[i] = newVolt;
666 }
667 if (busType[i] == BUS_PV) {
668 std::complex<double> yeSum(0.0, 0.0);
669 for (int k = 0; k < m_numberOfBuses; k++) {
670 if (i != k) {
671 // Sum { Y[i,k] * E[k] } | k = 1->n; k diff i
672 yeSum += m_yBus[i][k] * voltage[k];
673 }
674 }
675 std::complex<double> yeSumT = yeSum + (m_yBus[i][i] * voltage[i]);
676
677 // Q[i] = - Im( E*[i] * Sum { Y[i,k] * E[k] } )
678 std::complex<double> qCalc = std::conj(voltage[i]) * yeSumT;
679 power[i] = std::complex<double>(power[i].real(), -qCalc.imag());
680
681 // E[i] = (1/Y[i,i])*((P[i]-jQ[i])/E*[i] - Sum { Y[i,k] * E[k] (k diff i) })
682 std::complex<double> newVolt = (1.0 / m_yBus[i][i]) * (std::conj(power[i]) / std::conj(voltage[i]) - yeSum);
683
684 // Apply the acceleration factor.
685 newVolt = std::complex<double>(accFactor * (newVolt.real() - voltage[i].real()) + voltage[i].real(),
686 accFactor * (newVolt.imag() - voltage[i].imag()) + voltage[i].imag());
687
688 // Keep the same voltage magnitude.
689 voltage[i] = std::complex<double>(std::abs(voltage[i]) * std::cos(std::arg(newVolt)),
690 std::abs(voltage[i]) * std::sin(std::arg(newVolt)));
691 }
692
693 double busError = std::max(std::abs(voltage[i].real() - oldVoltage[i].real()),
694 std::abs(voltage[i].imag() - oldVoltage[i].imag()));
695
696 if (busError > error) error = busError;
697 }
698 return error;
699}
700
701void PowerFlow::NewtonRaphson(std::vector<BusType> busType,
702 std::vector<std::complex<double> >& voltage,
703 std::vector<std::complex<double> > power,
704 int numPV,
705 int numPQ,
706 std::vector<double> dPdQ,
707 double inertia)
708{
709 // Jacobian matrix
710 std::vector<std::vector<double> > jacobMatrix = CalculateJacobianMatrix(voltage, power, busType, numPV, numPQ);
711
712 // Apply inertia
713 for (unsigned int i = 0; i < dPdQ.size(); ++i) { dPdQ[i] = inertia * dPdQ[i]; }
714
715 // Calculate DeltaTheta DeltaV array
716 std::vector<double> dTdV = GaussianElimination(jacobMatrix, dPdQ);
717
718 // Update voltage array
719 int indexDT = 0;
720 int indexDV = numPQ + numPV;
721 for (int i = 0; i < m_numberOfBuses; i++) {
722 if (busType[i] != BUS_SLACK) {
723 if (busType[i] == BUS_PV) {
724 double newV = std::abs(voltage[i]);
725 double newT = std::arg(voltage[i]) + dTdV[indexDT];
726 voltage[i] = std::complex<double>(newV * std::cos(newT), newV * std::sin(newT));
727 indexDT++;
728 }
729 else {
730 //double newV = std::abs(voltage[i]) * (1.0 + dTdV[indexDV]);
731 double newV = std::abs(voltage[i]) + dTdV[indexDV];
732 double newT = std::arg(voltage[i]) + dTdV[indexDT];
733 voltage[i] = std::complex<double>(newV * std::cos(newT), newV * std::sin(newT));
734 indexDV++;
735 indexDT++;
736 }
737 }
738 }
739}
740
741std::vector<std::vector<double>>
742PowerFlow::CalculateJacobianMatrix(
743 const std::vector<std::complex<double>>& voltage,
744 const std::vector<std::complex<double>>& power,
745 const std::vector<BusType>& busType,
746 int numPV,
747 int numPQ)
748{
749 int nvar = 2 * numPQ + numPV;
750
751 std::vector<std::vector<double>> J(nvar, std::vector<double>(nvar, 0.0));
752
753 int n = m_numberOfBuses;
754
755 // Precompute magnitudes and angles
756 std::vector<double> V(n);
757 std::vector<double> T(n);
758
759 for (int i = 0; i < n; i++)
760 {
761 V[i] = std::abs(voltage[i]);
762 T[i] = std::arg(voltage[i]);
763 }
764
765 // Compute currents
766 std::vector<std::complex<double>> I(n);
767
768 for (int i = 0; i < n; i++)
769 {
770 I[i] = 0;
771 for (int k = 0; k < n; k++)
772 I[i] += m_yBus[i][k] * voltage[k];
773 }
774
775 // Compute power injections
776 std::vector<std::complex<double>> S(n);
777
778 for (int i = 0; i < n; i++)
779 S[i] = voltage[i] * std::conj(I[i]);
780
781 int rowP = 0;
782 int rowQ = numPV + numPQ;
783
784 for (int i = 0; i < n; i++)
785 {
786 if (busType[i] == BUS_SLACK)
787 continue;
788
789 int colTheta = 0;
790 int colV = numPV + numPQ;
791
792 for (int j = 0; j < n; j++)
793 {
794 if (busType[j] == BUS_SLACK)
795 continue;
796
797 double G = m_yBus[i][j].real();
798 double B = m_yBus[i][j].imag();
799
800 double dtheta = T[i] - T[j];
801
802 if (i == j)
803 {
804 double Pi = S[i].real();
805 double Qi = S[i].imag();
806
807 // J11
808 J[rowP][colTheta] = -Qi - B * V[i] * V[i];
809
810 // J12 only PQ
811 if (busType[i] == BUS_PQ)
812 J[rowQ][colTheta] = Pi - G * V[i] * V[i];
813 }
814 else
815 {
816 // J11
817 J[rowP][colTheta] = V[i] * V[j] * (G * sin(dtheta) - B * cos(dtheta));
818
819 if (busType[i] == BUS_PQ)// J21
820 J[rowQ][colTheta] = -V[i] * V[j] * (G * cos(dtheta) + B * sin(dtheta));
821 }
822
823 colTheta++;
824 }
825
826 colTheta = 0;
827
828 for (int j = 0; j < n; j++)
829 {
830 if (busType[j] != BUS_PQ)
831 continue;
832
833 double G = m_yBus[i][j].real();
834 double B = m_yBus[i][j].imag();
835
836 double dtheta = T[i] - T[j];
837
838 if (i == j)
839 {
840 double Pi = S[i].real();
841 double Qi = S[i].imag();
842
843 // J12
844 J[rowP][colV] = Pi / V[i] + G * V[i];
845
846 if (busType[i] == BUS_PQ) // J22
847 J[rowQ][colV] = Qi / V[i] - B * V[i];
848 }
849 else
850 {
851 // J12
852 J[rowP][colV] = V[i] * (G * cos(dtheta) + B * sin(dtheta));
853
854 if (busType[i] == BUS_PQ) // J22
855 J[rowQ][colV] = V[i] * (G * sin(dtheta) - B * cos(dtheta));
856 }
857 colV++;
858 }
859 rowP++;
860
861 if (busType[i] == BUS_PQ)
862 rowQ++;
863 }
864
865 return J;
866}
867
868bool PowerFlow::CheckReactiveLimits(std::vector<BusType>& busType,
869 std::vector<ReactiveLimits>& reactiveLimit,
870 std::vector<std::complex<double> >& power,
871 std::vector<std::complex<double> > loadPower)
872{
873 bool limitReach = false;
874 for (int i = 0; i < m_numberOfBuses; i++) {
875 if (busType[i] == BUS_PV) {
876 if (reactiveLimit[i].maxLimitType == RL_LIMITED) {
877 if (power[i].imag() - loadPower[i].imag() > reactiveLimit[i].maxLimit) {
878 power[i] = std::complex<double>(power[i].real(), reactiveLimit[i].maxLimit + loadPower[i].imag());
879 busType[i] = BUS_PQ;
880 reactiveLimit[i].limitReached = RL_MAX_REACHED;
881 limitReach = true;
882 }
883 }
884 if (reactiveLimit[i].minLimitType == RL_LIMITED) {
885 if (power[i].imag() - loadPower[i].imag() < reactiveLimit[i].minLimit) {
886 power[i] = std::complex<double>(power[i].real(), reactiveLimit[i].minLimit + loadPower[i].imag());
887 busType[i] = BUS_PQ;
888 reactiveLimit[i].limitReached = RL_MIN_REACHED;
889 limitReach = true;
890 }
891 }
892 }
893 }
894 return limitReach;
895}
896
897bool PowerFlow::CalculateMotorsReactivePower(std::vector<std::complex<double> > voltage,
898 std::vector<std::complex<double> >& power)
899{
900 for (unsigned int i = 0; i < m_indMotorList.size(); ++i) {
901 IndMotor* motor = m_indMotorList[i];
902 auto data = motor->GetElectricalData();
903 if (motor->IsOnline() && data.calcQInPowerFlow) {
904 double oldQ = data.qValue;
905 if (!motor->CalculateReactivePower(std::abs(voltage[data.busNum]))) return false;
906 double dQ = oldQ - motor->GetElectricalData().qValue;
907 power[data.busNum] += std::complex<double>(0.0, dQ);
908 }
909 }
910 return true;
911}
BUS_SLACK
@ BUS_SLACK
Definition ElectricCalculation.h:45
RL_MAX_REACHED
@ RL_MAX_REACHED
Definition ElectricCalculation.h:58
RL_MIN_REACHED
@ RL_MIN_REACHED
Definition ElectricCalculation.h:59
RL_LIMITED
@ RL_LIMITED
Definition ElectricCalculation.h:56
RL_UNLIMITED_SOURCE
@ RL_UNLIMITED_SOURCE
Definition ElectricCalculation.h:57
Bus
Node for power elements. All others power elements are connected through this.
Definition Bus.h:86
EMTElement
Element to connect ATP-EMTP.
Definition EMTElement.h:65
ElectricCalculation
Base class for electrical calculations providing general utility methods.
Definition ElectricCalculation.h:116
ElectricCalculation::m_busList
std::vector< Bus * > m_busList
List of buses in the system.
Definition ElectricCalculation.h:370
ElectricCalculation::UpdateElementsPowerFlow
virtual void UpdateElementsPowerFlow(std::vector< std::complex< double > > voltage, std::vector< std::complex< double > > power, std::vector< BusType > busType, std::vector< ReactiveLimits > reactiveLimit, double systemPowerBase)
Update the elements after the power flow calculation.
Definition ElectricCalculation.cpp:453
ElectricCalculation::m_loadList
std::vector< Load * > m_loadList
List of load elements in the system.
Definition ElectricCalculation.h:380
ElectricCalculation::CalculateEMTElementsAdmittance
bool CalculateEMTElementsAdmittance(const double &basePower, wxString &errorMsg)
Calculate the admittance of EMT elements.
Definition ElectricCalculation.cpp:1057
ElectricCalculation::CalculateEMTPowerError
double CalculateEMTPowerError(const std::vector< std::complex< double > > &voltage, std::vector< std::complex< double > > &power, const double &basePower, wxString &errorMsg)
Calculate the power mismatch error for EMT simulation.
Definition ElectricCalculation.cpp:1109
ElectricCalculation::m_indMotorList
std::vector< IndMotor * > m_indMotorList
List of induction motors in the system.
Definition ElectricCalculation.h:374
ElectricCalculation::CalculateEMTElementsPower
bool CalculateEMTElementsPower(const double &basePower, wxString &errorMsg, bool updateCurrent=true)
Calculate the power of EMT elements.
Definition ElectricCalculation.cpp:1075
ElectricCalculation::m_syncGeneratorList
std::vector< SyncGenerator * > m_syncGeneratorList
List of synchronous generators in the system.
Definition ElectricCalculation.h:382
ElectricCalculation::m_emtElementList
std::vector< EMTElement * > m_emtElementList
List of electromagnetic transient (EMT) elements in the system.
Definition ElectricCalculation.h:390
ElectricCalculation::GaussianElimination
std::vector< std::complex< double > > GaussianElimination(std::vector< std::vector< std::complex< double > > > matrix, std::vector< std::complex< double > > array)
Solve a linear system using Gaussian elimination (complex version).
Definition ElectricCalculation.cpp:842
ElectricCalculation::GetYBus
virtual bool GetYBus(std::vector< std::vector< std::complex< double > > > &yBus, double systemPowerBase, YBusSequence sequence=POSITIVE_SEQ, bool includeSyncMachines=false, bool allLoadsAsImpedances=false, bool usePowerFlowVoltagesOnImpedances=false)
Get the admittance matrix from the list of elements (use GetElementsFromList first).
Definition ElectricCalculation.cpp:81
ElectricCalculation::m_syncMotorList
std::vector< SyncMotor * > m_syncMotorList
List of synchronous motors in the system.
Definition ElectricCalculation.h:384
Element::GetParentList
virtual std::vector< Element * > GetParentList() const
Get the parent list.
Definition Element.h:559
Element::IsOnline
bool IsOnline() const
Checks if the element is online or offline.
Definition Element.h:226
IndMotor
Induction motor power element.
Definition IndMotor.h:119
Load
Loas shunt power element.
Definition Load.h:74
SyncGenerator
Synchronous generator power element.
Definition SyncGenerator.h:146
SyncMotor
Synchronous motor (synchronous compensator) power element.
Definition SyncMotor.h:135