Power System Platform  2026w34a-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 m_tapAdjustmentsCount = 0;
247 return true;
248}
249
250bool PowerFlow::RunGaussSeidel(double systemPowerBase,
251 int maxIteration,
252 double error,
253 double initAngle,
254 double accFactor)
255{
256 std::vector<BusType> busType; // Bus type
257 std::vector<std::complex<double> > voltage; // Voltage of buses
258 std::vector<std::complex<double> > power; // Injected power
259 std::vector<std::complex<double> > loadPower; // Only the load power
260 std::vector<ReactiveLimits> reactiveLimit; // Limit of reactive power on PV buses
261
262 if (!InitPowerFlow(busType, voltage, power, loadPower, reactiveLimit, systemPowerBase, initAngle)) return false;
263
264 // Gauss-Seidel method
265 std::vector<std::complex<double> > oldVoltage; // Old voltage array.
266 oldVoltage.resize(voltage.size());
267
268 auto oldBusType = busType;
269
270 int iteration = 0; // Current itaration number.
271 double emtPowerError = 1e3; // EMT Elements power error.
272
273 while (true) {
274 // Reach the max number of iterations.
275 if (iteration >= maxIteration) {
276 m_errorMsg = _("The maximum number of iterations was reached.");
277 return false;
278 }
279
280 // Calculate induction motor reactive power
281 if (!CalculateMotorsReactivePower(voltage, power)) {
282 m_errorMsg = _("It was not possible to solve the induction motors.");
283 return false;
284 }
285
286 // Update the old voltage array to current iteration values.
287 for (int i = 0; i < m_numberOfBuses; i++) oldVoltage[i] = voltage[i];
288
289 double iterationError = GaussSeidel(busType, voltage, oldVoltage, power, accFactor);
290
291 if (HasInvalidValue(voltage)) {
292 m_errorMsg = _("The power flow solution has invalid voltage values.");
293 ResetVoltages();
294 return false;
295 }
296 if (HasInvalidValue(power)) {
297 m_errorMsg = _("The power flow solution has invalid power values.");
298 ResetVoltages();
299 return false;
300 }
301
302 if (iterationError < error) {
303 // Calculate power error in EMT elements.
304 if (!m_emtElementList.empty()) {
305 double oldError = emtPowerError;
306 emtPowerError = CalculateEMTPowerError(voltage, power, systemPowerBase, m_errorMsg);
307 //wxMessageBox(wxString::Format(wxT("It %d, EMT Power Error: %e (base = %e)"), iteration, emtPowerError, error));
308 if (abs(emtPowerError - oldError) < 1e-12) { // Without change in error the convergence is impossible
309 m_errorMsg = _("Impossible to reach a convergence with the current Electromagnetic Transient Elements.");
310 return false;
311 }
312 }
313 else emtPowerError = 0.0;
314
315 bool qLimitReached = CheckReactiveLimits(busType, reactiveLimit, power, loadPower);
316 bool tapAdjusted = AdjustTapChangers(voltage, systemPowerBase);
317
318 if (!qLimitReached && !tapAdjusted && emtPowerError < error) break;
319 }
320
321 iteration++;
322 }
323 m_iterations = iteration;
324
325 // Adjust the power array.
326 for (int i = 0; i < m_numberOfBuses; i++) {
327 std::complex<double> sBus = std::complex<double>(0.0, 0.0);
328 for (int j = 0; j < m_numberOfBuses; j++) sBus += voltage[i] * std::conj(voltage[j]) * std::conj(m_yBus[i][j]);
329 power[i] = sBus;
330 }
331
332
333 UpdateElementsPowerFlow(voltage, power, oldBusType, reactiveLimit, systemPowerBase);
334
335 return true;
336}
337
338bool PowerFlow::HasInvalidValue(const std::vector< std::complex<double> >& value)
339{
340 for (size_t i = 0; i < value.size(); i++)
341 {
342 if (!std::isfinite(value[i].real()) ||
343 !std::isfinite(value[i].imag()))
344 {
345 return true;
346 }
347 }
348 return false;
349}
350
351bool PowerFlow::RunNewtonRaphson(double systemPowerBase,
352 int maxIteration,
353 double error,
354 double initAngle,
355 double inertia)
356{
357 std::vector<BusType> busType; // Bus type
358 std::vector<std::complex<double> > voltage; // Voltage of buses
359 std::vector<std::complex<double> > power; // Injected power
360 std::vector<std::complex<double> > loadPower; // Only the load power
361 std::vector<ReactiveLimits> reactiveLimit; // Limit of reactive power on PV buses
362
363 if (!InitPowerFlow(busType, voltage, power, loadPower, reactiveLimit, systemPowerBase, initAngle)) return false;
364 auto oldBusType = busType;
365
366 // Newton-Raphson method
367 int numPQ = 0; // Number of PQ buses
368 int numPV = 0; // Number of PV buses
369 GetNumPVPQ(busType, numPQ, numPV);
370
371 // DeltaP and DeltaQ array
372 std::vector<double> dPdQ;
373 dPdQ.resize(numPV + 2 * numPQ, 0.0);
374
375 int iteration = 0; // Current iteration number.
376 while (true) {
377 // Reach the max number of iterations.
378 if (iteration >= maxIteration) {
379 m_errorMsg = _("The maximum number of iterations was reached.");
380 return false;
381 }
382
383 // Calculate induction motor reactive power
384 if (!CalculateMotorsReactivePower(voltage, power)) {
385 m_errorMsg = _("It was not possible to solve the induction motors.");
386 return false;
387 }
388
389 // Calculate dPdQ array
390
391 // Fill it with zeros
392 std::fill(dPdQ.begin(), dPdQ.end(), 0.0);
393
394 int indexDP = 0;
395 int indexDQ = numPQ + numPV;
396 for (int i = 0; i < m_numberOfBuses; i++) {
397 if (busType[i] != BUS_SLACK) {
398 for (int j = 0; j < m_numberOfBuses; j++) {
399 // PV ou PQ bus
400 std::complex<double> sInj = std::conj(m_yBus[i][j]) * voltage[i] * std::conj(voltage[j]);
401 dPdQ[indexDP] += sInj.real();
402
403 // PQ bus
404 if (busType[i] == BUS_PQ) dPdQ[indexDQ] += sInj.imag();
405 }
406
407 // PQ or PV bus
408 dPdQ[indexDP] = power[i].real() - dPdQ[indexDP];
409 indexDP++;
410
411 // PQ bus
412 if (busType[i] == BUS_PQ) {
413 dPdQ[indexDQ] = power[i].imag() - dPdQ[indexDQ];
414 indexDQ++;
415 }
416 }
417 }
418
419 // Calculate the iteration error
420 double iterationError = 0.0;
421 for (unsigned int i = 0; i < dPdQ.size(); ++i) {
422 if (iterationError < std::abs(dPdQ[i])) iterationError = std::abs(dPdQ[i]);
423 }
424
425 // Check if the iteration error is less than tolerance, also check if any reactive limit was reached.
426 // If any reactive limit was reached, change the bus type.
427 if (iterationError < error) {
428 for (int i = 0; i < m_numberOfBuses; i++) {
429 std::complex<double> sBus(0.0, 0.0);
430 for (int j = 0; j < m_numberOfBuses; j++)
431 sBus += voltage[i] * std::conj(voltage[j]) * std::conj(m_yBus[i][j]);
432 power[i] = sBus;
433 }
434
435 double emtPowerError = CalculateEMTPowerError(voltage, power, systemPowerBase, m_errorMsg);
436
437 bool qLimitReached = CheckReactiveLimits(busType, reactiveLimit, power, loadPower);
438 bool tapAdjusted = AdjustTapChangers(voltage, systemPowerBase);
439
440 if (!qLimitReached && !tapAdjusted && emtPowerError < error)
441 break;
442 else {
443 GetNumPVPQ(busType, numPQ, numPV);
444 dPdQ.assign(numPV + 2 * numPQ, 0.0);
445 }
446 }
447
448 NewtonRaphson(busType, voltage, power, numPV, numPQ, dPdQ, inertia);
449
450 if (HasInvalidValue(voltage)) {
451 m_errorMsg = _("The power flow solution has invalid voltage values.");
452 ResetVoltages();
453 return false;
454 }
455 if (HasInvalidValue(power)) {
456 m_errorMsg = _("The power flow solution has invalid power values.");
457 ResetVoltages();
458 return false;
459 }
460
461 iteration++;
462 }
463 m_iterations = iteration;
464
465
466
467 // Adjust the power array.
468 for (int i = 0; i < m_numberOfBuses; i++) {
469 std::complex<double> sBus = std::complex<double>(0.0, 0.0);
470 for (int j = 0; j < m_numberOfBuses; j++) sBus += voltage[i] * std::conj(voltage[j]) * std::conj(m_yBus[i][j]);
471 power[i] = sBus;
472 }
473
474 UpdateElementsPowerFlow(voltage, power, oldBusType, reactiveLimit, systemPowerBase);
475
476 return true;
477}
478
479bool PowerFlow::RunGaussNewton(double systemPowerBase,
480 int maxIteration,
481 double error,
482 double initAngle,
483 double accFactor,
484 double gaussTol,
485 double inertia)
486{
487 std::vector<BusType> busType; // Bus type
488 std::vector<std::complex<double> > voltage; // Voltage of buses
489 std::vector<std::complex<double> > power; // Injected power
490 std::vector<std::complex<double> > loadPower; // Only the load power
491 std::vector<ReactiveLimits> reactiveLimit; // Limit of reactive power on PV buses
492
493 if (!InitPowerFlow(busType, voltage, power, loadPower, reactiveLimit, systemPowerBase, initAngle)) return false;
494
495 // Gauss-Seidel method
496 std::vector<std::complex<double> > oldVoltage; // Old voltage array.
497 oldVoltage.resize(voltage.size());
498
499 auto oldBusType = busType;
500
501 int iteration = 0; // Current itaration number.
502
503 while (true) {
504 // Reach the max number of iterations.
505 if (iteration >= maxIteration) {
506 m_errorMsg = _("The maximum number of iterations was reached.");
507 return false;
508 }
509
510 // Calculate induction motor reactive power
511 if (!CalculateMotorsReactivePower(voltage, power)) {
512 m_errorMsg = _("It was not possible to solve the induction motors.");
513 return false;
514 }
515
516 // Update the old voltage array to current iteration values.
517 for (int i = 0; i < m_numberOfBuses; i++) oldVoltage[i] = voltage[i];
518
519 double iterationError = GaussSeidel(busType, voltage, oldVoltage, power, accFactor);
520
521 if (iterationError < gaussTol) break;
522
523 iteration++;
524 }
525
526 // Newton-Raphson method
527 int numPQ = 0; // Number of PQ buses
528 int numPV = 0; // Number of PV buses
529 GetNumPVPQ(busType, numPQ, numPV);
530
531 // DeltaP and DeltaQ array
532 std::vector<double> dPdQ;
533 dPdQ.resize(numPV + 2 * numPQ, 0.0);
534
535 while (true) {
536 // Reach the max number of iterations.
537 if (iteration >= maxIteration) {
538 m_errorMsg = _("The maximum number of iterations was reached.");
539 return false;
540 }
541
542 // Calculate induction motor reactive power
543 if (!CalculateMotorsReactivePower(voltage, power)) {
544 m_errorMsg = _("It was not possible to solve the induction motors.");
545 return false;
546 }
547
548 // Calculate dPdQ array
549
550 // Fill it with zeros
551 std::fill(dPdQ.begin(), dPdQ.end(), 0.0);
552
553 int indexDP = 0;
554 int indexDQ = numPQ + numPV;
555 for (int i = 0; i < m_numberOfBuses; i++) {
556 if (busType[i] != BUS_SLACK) {
557 for (int j = 0; j < m_numberOfBuses; j++) {
558 // PV ou PQ bus
559 std::complex<double> sInj = std::conj(m_yBus[i][j]) * voltage[i] * std::conj(voltage[j]);
560 dPdQ[indexDP] += sInj.real();
561
562 // PQ bus
563 if (busType[i] == BUS_PQ) dPdQ[indexDQ] += sInj.imag();
564 }
565
566 // PQ or PV bus
567 dPdQ[indexDP] = power[i].real() - dPdQ[indexDP];
568 indexDP++;
569
570 // PQ bus
571 if (busType[i] == BUS_PQ) {
572 dPdQ[indexDQ] = power[i].imag() - dPdQ[indexDQ];
573 indexDQ++;
574 }
575 }
576 }
577
578 // Calculate the iteration error
579 double iterationError = 0.0;
580 for (unsigned int i = 0; i < dPdQ.size(); ++i) {
581 if (iterationError < std::abs(dPdQ[i])) iterationError = std::abs(dPdQ[i]);
582 }
583
584 // Check if the iteration error is less than tolerance, also check if any reactive limit was reached.
585 // If any reactive limit was reached, change the bus type.
586 if (iterationError < error) {
587 double emtPowerError = CalculateEMTPowerError(voltage, power, systemPowerBase, m_errorMsg);
588
589 bool qLimitReached = CheckReactiveLimits(busType, reactiveLimit, power, loadPower);
590 bool tapAdjusted = AdjustTapChangers(voltage, systemPowerBase);
591
592 if (!qLimitReached && !tapAdjusted && emtPowerError < error)
593 break;
594 else {
595 GetNumPVPQ(busType, numPQ, numPV);
596 dPdQ.clear();
597 dPdQ.resize(numPV + 2 * numPQ, 0.0);
598 }
599 }
600
601 NewtonRaphson(busType, voltage, power, numPV, numPQ, dPdQ, inertia);
602
603 iteration++;
604 }
605 m_iterations = iteration;
606
607 // Adjust the power array.
608 for (int i = 0; i < m_numberOfBuses; i++) {
609 std::complex<double> sBus = std::complex<double>(0.0, 0.0);
610 for (int j = 0; j < m_numberOfBuses; j++) sBus += voltage[i] * std::conj(voltage[j]) * std::conj(m_yBus[i][j]);
611 power[i] = sBus;
612 }
613
614 UpdateElementsPowerFlow(voltage, power, oldBusType, reactiveLimit, systemPowerBase);
615
616 return true;
617}
618
619void PowerFlow::ResetVoltages()
620{
621 for (auto* bus : m_busList) {
622 BusElectricalData data = bus->GetElectricalData();
623 if (bus->IsOnline()) {
624 data.voltage = std::complex<double>(1.0, 0.0);
625 if (data.isVoltageControlled) {
626 data.voltage = std::complex<double>(data.controlledVoltage, 0.0);
627 }
628 }
629 else
630 data.voltage = std::complex<double>(0.0, 0.0);
631 data.harmonicOrder.clear();
632 data.harmonicVoltage.clear();
633 data.thd = 0.0;
634 bus->SetElectricalData(data);
635 }
636 for (auto* transf : m_transformerList) {
637 auto data = transf->GetElectricalData();
638 if (data.hasTapChanger && data.nominalTurnsRatio > 0.0) {
639 data.turnsRatio = data.nominalTurnsRatio;
640 transf->SetElectricaData(data);
641 }
642 }
643}
644
645void PowerFlow::GetNumPVPQ(std::vector<BusType> busType, int& numPQ, int& numPV)
646{
647 numPQ = 0;
648 numPV = 0;
649 for (auto it = busType.begin(), itEnd = busType.end(); it != itEnd; ++it) {
650 if (*it == BUS_PQ)
651 numPQ++;
652 else if (*it == BUS_PV)
653 numPV++;
654 }
655}
656
657double PowerFlow::GaussSeidel(std::vector<BusType> busType,
658 std::vector<std::complex<double> >& voltage,
659 std::vector<std::complex<double> > oldVoltage,
660 std::vector<std::complex<double> >& power,
661 double accFactor)
662{
663 double error = 0.0;
664
665 for (int i = 0; i < m_numberOfBuses; i++) {
666 if (busType[i] == BUS_PQ) {
667 std::complex<double> yeSum(0.0, 0.0);
668 for (int k = 0; k < m_numberOfBuses; k++) {
669 if (i != k) {
670 // Sum { Y[i,k] * E[k] } | k = 1->n; k diff i
671 yeSum += m_yBus[i][k] * voltage[k];
672 }
673 }
674
675 // E[i] = (1/Y[i,i])*((P[i]-jQ[i])/E*[i] - Sum { Y[i,k] * E[k] (k diff i) })
676 std::complex<double> newVolt = (1.0 / m_yBus[i][i]) * (std::conj(power[i]) / std::conj(voltage[i]) - yeSum);
677
678 // Apply the acceleration factor.
679 newVolt = std::complex<double>(accFactor * (newVolt.real() - voltage[i].real()) + voltage[i].real(),
680 accFactor * (newVolt.imag() - voltage[i].imag()) + voltage[i].imag());
681
682 voltage[i] = newVolt;
683 }
684 if (busType[i] == BUS_PV) {
685 std::complex<double> yeSum(0.0, 0.0);
686 for (int k = 0; k < m_numberOfBuses; k++) {
687 if (i != k) {
688 // Sum { Y[i,k] * E[k] } | k = 1->n; k diff i
689 yeSum += m_yBus[i][k] * voltage[k];
690 }
691 }
692 std::complex<double> yeSumT = yeSum + (m_yBus[i][i] * voltage[i]);
693
694 // Q[i] = - Im( E*[i] * Sum { Y[i,k] * E[k] } )
695 std::complex<double> qCalc = std::conj(voltage[i]) * yeSumT;
696 power[i] = std::complex<double>(power[i].real(), -qCalc.imag());
697
698 // E[i] = (1/Y[i,i])*((P[i]-jQ[i])/E*[i] - Sum { Y[i,k] * E[k] (k diff i) })
699 std::complex<double> newVolt = (1.0 / m_yBus[i][i]) * (std::conj(power[i]) / std::conj(voltage[i]) - yeSum);
700
701 // Apply the acceleration factor.
702 newVolt = std::complex<double>(accFactor * (newVolt.real() - voltage[i].real()) + voltage[i].real(),
703 accFactor * (newVolt.imag() - voltage[i].imag()) + voltage[i].imag());
704
705 // Keep the same voltage magnitude.
706 voltage[i] = std::complex<double>(std::abs(voltage[i]) * std::cos(std::arg(newVolt)),
707 std::abs(voltage[i]) * std::sin(std::arg(newVolt)));
708 }
709
710 double busError = std::max(std::abs(voltage[i].real() - oldVoltage[i].real()),
711 std::abs(voltage[i].imag() - oldVoltage[i].imag()));
712
713 if (busError > error) error = busError;
714 }
715 return error;
716}
717
718void PowerFlow::NewtonRaphson(std::vector<BusType> busType,
719 std::vector<std::complex<double> >& voltage,
720 std::vector<std::complex<double> > power,
721 int numPV,
722 int numPQ,
723 std::vector<double> dPdQ,
724 double inertia)
725{
726 // Jacobian matrix
727 std::vector<std::vector<double> > jacobMatrix = CalculateJacobianMatrix(voltage, power, busType, numPV, numPQ);
728
729 // Apply inertia
730 for (unsigned int i = 0; i < dPdQ.size(); ++i) { dPdQ[i] = inertia * dPdQ[i]; }
731
732 // Calculate DeltaTheta DeltaV array
733 std::vector<double> dTdV = GaussianElimination(jacobMatrix, dPdQ);
734
735 // Update voltage array
736 int indexDT = 0;
737 int indexDV = numPQ + numPV;
738 for (int i = 0; i < m_numberOfBuses; i++) {
739 if (busType[i] != BUS_SLACK) {
740 if (busType[i] == BUS_PV) {
741 double newV = std::abs(voltage[i]);
742 double newT = std::arg(voltage[i]) + dTdV[indexDT];
743 voltage[i] = std::complex<double>(newV * std::cos(newT), newV * std::sin(newT));
744 indexDT++;
745 }
746 else {
747 //double newV = std::abs(voltage[i]) * (1.0 + dTdV[indexDV]);
748 double newV = std::abs(voltage[i]) + dTdV[indexDV];
749 double newT = std::arg(voltage[i]) + dTdV[indexDT];
750 voltage[i] = std::complex<double>(newV * std::cos(newT), newV * std::sin(newT));
751 indexDV++;
752 indexDT++;
753 }
754 }
755 }
756}
757
758std::vector<std::vector<double>>
759PowerFlow::CalculateJacobianMatrix(
760 const std::vector<std::complex<double>>& voltage,
761 const std::vector<std::complex<double>>& power,
762 const std::vector<BusType>& busType,
763 int numPV,
764 int numPQ)
765{
766 int nvar = 2 * numPQ + numPV;
767
768 std::vector<std::vector<double>> J(nvar, std::vector<double>(nvar, 0.0));
769
770 int n = m_numberOfBuses;
771
772 // Precompute magnitudes and angles
773 std::vector<double> V(n);
774 std::vector<double> T(n);
775
776 for (int i = 0; i < n; i++)
777 {
778 V[i] = std::abs(voltage[i]);
779 T[i] = std::arg(voltage[i]);
780 }
781
782 // Compute currents
783 std::vector<std::complex<double>> I(n);
784
785 for (int i = 0; i < n; i++)
786 {
787 I[i] = 0;
788 for (int k = 0; k < n; k++)
789 I[i] += m_yBus[i][k] * voltage[k];
790 }
791
792 // Compute power injections
793 std::vector<std::complex<double>> S(n);
794
795 for (int i = 0; i < n; i++)
796 S[i] = voltage[i] * std::conj(I[i]);
797
798 int rowP = 0;
799 int rowQ = numPV + numPQ;
800
801 for (int i = 0; i < n; i++)
802 {
803 if (busType[i] == BUS_SLACK)
804 continue;
805
806 int colTheta = 0;
807 int colV = numPV + numPQ;
808
809 for (int j = 0; j < n; j++)
810 {
811 if (busType[j] == BUS_SLACK)
812 continue;
813
814 double G = m_yBus[i][j].real();
815 double B = m_yBus[i][j].imag();
816
817 double dtheta = T[i] - T[j];
818
819 if (i == j)
820 {
821 double Pi = S[i].real();
822 double Qi = S[i].imag();
823
824 // J11
825 J[rowP][colTheta] = -Qi - B * V[i] * V[i];
826
827 // J12 only PQ
828 if (busType[i] == BUS_PQ)
829 J[rowQ][colTheta] = Pi - G * V[i] * V[i];
830 }
831 else
832 {
833 // J11
834 J[rowP][colTheta] = V[i] * V[j] * (G * sin(dtheta) - B * cos(dtheta));
835
836 if (busType[i] == BUS_PQ)// J21
837 J[rowQ][colTheta] = -V[i] * V[j] * (G * cos(dtheta) + B * sin(dtheta));
838 }
839
840 colTheta++;
841 }
842
843 colTheta = 0;
844
845 for (int j = 0; j < n; j++)
846 {
847 if (busType[j] != BUS_PQ)
848 continue;
849
850 double G = m_yBus[i][j].real();
851 double B = m_yBus[i][j].imag();
852
853 double dtheta = T[i] - T[j];
854
855 if (i == j)
856 {
857 double Pi = S[i].real();
858 double Qi = S[i].imag();
859
860 // J12
861 J[rowP][colV] = Pi / V[i] + G * V[i];
862
863 if (busType[i] == BUS_PQ) // J22
864 J[rowQ][colV] = Qi / V[i] - B * V[i];
865 }
866 else
867 {
868 // J12
869 J[rowP][colV] = V[i] * (G * cos(dtheta) + B * sin(dtheta));
870
871 if (busType[i] == BUS_PQ) // J22
872 J[rowQ][colV] = V[i] * (G * sin(dtheta) - B * cos(dtheta));
873 }
874 colV++;
875 }
876 rowP++;
877
878 if (busType[i] == BUS_PQ)
879 rowQ++;
880 }
881
882 return J;
883}
884
885bool PowerFlow::CheckReactiveLimits(std::vector<BusType>& busType,
886 std::vector<ReactiveLimits>& reactiveLimit,
887 std::vector<std::complex<double> >& power,
888 std::vector<std::complex<double> > loadPower)
889{
890 bool limitReach = false;
891 for (int i = 0; i < m_numberOfBuses; i++) {
892 if (busType[i] == BUS_PV) {
893 if (reactiveLimit[i].maxLimitType == RL_LIMITED) {
894 if (power[i].imag() - loadPower[i].imag() > reactiveLimit[i].maxLimit) {
895 power[i] = std::complex<double>(power[i].real(), reactiveLimit[i].maxLimit + loadPower[i].imag());
896 busType[i] = BUS_PQ;
897 reactiveLimit[i].limitReached = RL_MAX_REACHED;
898 limitReach = true;
899 }
900 }
901 if (reactiveLimit[i].minLimitType == RL_LIMITED) {
902 if (power[i].imag() - loadPower[i].imag() < reactiveLimit[i].minLimit) {
903 power[i] = std::complex<double>(power[i].real(), reactiveLimit[i].minLimit + loadPower[i].imag());
904 busType[i] = BUS_PQ;
905 reactiveLimit[i].limitReached = RL_MIN_REACHED;
906 limitReach = true;
907 }
908 }
909 }
910 }
911 return limitReach;
912}
913
914bool PowerFlow::CalculateMotorsReactivePower(std::vector<std::complex<double> > voltage,
915 std::vector<std::complex<double> >& power)
916{
917 for (unsigned int i = 0; i < m_indMotorList.size(); ++i) {
918 IndMotor* motor = m_indMotorList[i];
919 auto data = motor->GetElectricalData();
920 if (motor->IsOnline() && data.calcQInPowerFlow) {
921 double oldQ = data.qValue;
922 if (!motor->CalculateReactivePower(std::abs(voltage[data.busNum]))) return false;
923 double dQ = oldQ - motor->GetElectricalData().qValue;
924 power[data.busNum] += std::complex<double>(0.0, dQ);
925 }
926 }
927 return true;
928}
929
930bool PowerFlow::AdjustTapChangers(const std::vector<std::complex<double> >& voltage, double systemPowerBase)
931{
932 if (m_tapAdjustmentsCount >= 40) {
933 // Limit total adjustments to avoid infinite cycling
934 return false;
935 }
936
937 bool anyAdjusted = false;
938
939 for (Transformer* transformer : m_transformerList) {
940 if (!transformer->IsOnline() || transformer->GetParentList().size() < 2) continue;
941
942 TransformerElectricalData data = transformer->GetElectricalData();
943 if (!data.hasTapChanger) continue;
944
945 Bus* bus1 = static_cast<Bus*>(transformer->GetParentList()[0]);
946 Bus* bus2 = static_cast<Bus*>(transformer->GetParentList()[1]);
947 if (!bus1 || !bus2) continue;
948
949 int n1 = bus1->GetElectricalData().number;
950 int n2 = bus2->GetElectricalData().number;
951
952 int nCtrl = (data.oltcControlledBus == 0) ? n1 : n2;
953 bool isPrimary = (data.oltcControlledBus == 0);
954
955 if (nCtrl < 0 || nCtrl >= (int)voltage.size()) continue;
956
957 double vCtrl = std::abs(voltage[nCtrl]);
958 double vTarget = data.oltcTargetVoltage;
959 double deadband = std::max(data.oltcVoltageDeadband, 1e-4);
960
961 if (data.oltcIsDiscrete && data.oltcTapStep > 1e-4) {
962 deadband = std::max(deadband, 0.5 * data.oltcTapStep * vCtrl);
963 }
964
965 double vDiff = vCtrl - vTarget;
966 if (std::abs(vDiff) <= deadband) continue;
967
968 double currentTap = data.turnsRatio;
969 if (currentTap <= 1e-4) currentTap = 1.0;
970
971
972 //tex:
973 // Model in PSP-UFU has ideal turns ratio $a$ on side 1 (Primary):
974 //
975 // $$\frac{V_1}{a} \approx V_2$$
976 // $$V_2 \approx \frac{V_1}{a}, \qquad V_1 \approx aV_2$$
977 //
978 // Secondary (side 2):
979 //
980 // $$\frac{dV_2}{da} \approx -\frac{V_2}{a}$$
981 //
982 // $$\Delta V_2 \approx -\frac{V_2}{a}\Delta a$$
983 // $$\Rightarrow\quad
984 // \Delta a \approx -\frac{a}{V_2}\Delta V_2
985 // = \frac{a}{V_2}(V_2 - V_{\mathrm{target}})
986 // = \frac{a}{V_2}v_{\mathrm{Diff}}$$
987 //
988 // Primary (side 1):
989 //
990 // $$\frac{dV_1}{da} \approx \frac{V_1}{a}$$
991 //
992 // $$\Delta V_1 \approx \frac{V_1}{a}\Delta a$$
993 // $$\Rightarrow\quad
994 // \Delta a \approx \frac{a}{V_1}\Delta V_1
995 // = -\frac{a}{V_1}(V_1 - V_{\mathrm{target}})
996 // = -\frac{a}{V_1}v_{\mathrm{Diff}}$$
997
998 double damping = 0.8;
999 double deltaTap = (isPrimary ? -1.0 : 1.0) * (currentTap / std::max(vCtrl, 0.1)) * vDiff * damping;
1000
1001 if (data.oltcIsDiscrete && data.oltcTapStep > 1e-4) {
1002 double step = data.oltcTapStep;
1003 double numSteps = std::round(deltaTap / step);
1004 if (numSteps == 0.0) {
1005 numSteps = (deltaTap > 0.0) ? 1.0 : -1.0;
1006 }
1007 deltaTap = numSteps * step;
1008 }
1009
1010 double newTap = currentTap + deltaTap;
1011
1012 // Limit to min/max range
1013 if (newTap < data.oltcMinTap) newTap = data.oltcMinTap;
1014 if (newTap > data.oltcMaxTap) newTap = data.oltcMaxTap;
1015
1016 if (data.oltcIsDiscrete && data.oltcTapStep > 1e-4) {
1017 double step = data.oltcTapStep;
1018 newTap = 1.0 + std::round((newTap - 1.0) / step) * step;
1019 if (newTap < data.oltcMinTap) newTap = data.oltcMinTap;
1020 if (newTap > data.oltcMaxTap) newTap = data.oltcMaxTap;
1021 }
1022
1023 if (std::abs(newTap - currentTap) > 1e-5) {
1024 data.turnsRatio = newTap;
1025 transformer->SetElectricaData(data);
1026 anyAdjusted = true;
1027 }
1028 }
1029
1030 if (anyAdjusted) {
1031 m_tapAdjustmentsCount++;
1032 GetYBus(m_yBus, systemPowerBase);
1033 }
1034
1035 return anyAdjusted;
1036}
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:87
EMTElement
Element to connect ATP-EMTP.
Definition EMTElement.h:83
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:1089
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:1141
ElectricCalculation::m_indMotorList
std::vector< IndMotor * > m_indMotorList
List of induction motors in the system.
Definition ElectricCalculation.h:374
ElectricCalculation::m_transformerList
std::vector< Transformer * > m_transformerList
List of transformers in the system.
Definition ElectricCalculation.h:386
ElectricCalculation::CalculateEMTElementsPower
bool CalculateEMTElementsPower(const double &basePower, wxString &errorMsg, bool updateCurrent=true)
Calculate the power of EMT elements.
Definition ElectricCalculation.cpp:1107
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:874
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:567
Element::IsOnline
bool IsOnline() const
Checks if the element is online or offline.
Definition Element.h:229
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
Transformer
Two-winding transformer power element with OLTC support.
Definition Transformer.h:96