Merge remote-tracking branch 'refs/remotes/origin/wip/currentState'

1 files changed, 427 insertions, 0 deletions

diff --git a/Project/ElectricCalculation.cpp b/Project/ElectricCalculation.cpp
new file mode 100644
index 0000000..3c4c1e5
--- /dev/null
+++ b/Project/ElectricCalculation.cpp
@@ -0,0 +1,427 @@
+#include "ElectricCalculation.h"
+
+ElectricCalculation::ElectricCalculation() {}
+ElectricCalculation::~ElectricCalculation() {}
+void ElectricCalculation::GetElementsFromList(std::vector<Element*> elementList)
+{
+    m_busList.clear();
+    m_capacitorList.clear();
+    m_indMotorList.clear();
+    m_inductorList.clear();
+    m_lineList.clear();
+    m_loadList.clear();
+    m_syncGeneratorList.clear();
+    m_syncMotorList.clear();
+    m_transformerList.clear();
+    // TODO: Bad design?
+    for(auto it = elementList.begin(); it != elementList.end(); it++) {
+        Element* element = *it;
+        if(typeid(*element) == typeid(Bus))
+            m_busList.push_back((Bus*)element);
+        else if(typeid(*element) == typeid(Capacitor))
+            m_capacitorList.push_back((Capacitor*)element);
+        else if(typeid(*element) == typeid(IndMotor))
+            m_indMotorList.push_back((IndMotor*)element);
+        else if(typeid(*element) == typeid(Inductor))
+            m_inductorList.push_back((Inductor*)element);
+        else if(typeid(*element) == typeid(Line))
+            m_lineList.push_back((Line*)element);
+        else if(typeid(*element) == typeid(Load))
+            m_loadList.push_back((Load*)element);
+        else if(typeid(*element) == typeid(SyncGenerator))
+            m_syncGeneratorList.push_back((SyncGenerator*)element);
+        else if(typeid(*element) == typeid(SyncMotor))
+            m_syncMotorList.push_back((SyncMotor*)element);
+        else if(typeid(*element) == typeid(Transformer))
+            m_transformerList.push_back((Transformer*)element);
+    }
+}
+
+bool ElectricCalculation::GetYBus(std::vector<std::vector<std::complex<double> > >& yBus, double systemPowerBase)
+{
+    if(m_busList.size() == 0) return false;
+
+    // Clear and fill with zeros the Ybus
+    yBus.clear();
+    for(int i = 0; i < (int)m_busList.size(); i++) {
+        std::vector<std::complex<double> > line;
+        for(int j = 0; j < (int)m_busList.size(); j++) {
+            line.push_back(std::complex<double>(0.0, 0.0));
+        }
+        yBus.push_back(line);
+    }
+
+    // Set buses numbers
+    int busNumber = 0;
+    for(auto itb = m_busList.begin(); itb != m_busList.end(); itb++) {
+        Bus* bus = *itb;
+        BusElectricalData data = bus->GetEletricalData();
+        data.number = busNumber;
+        bus->SetElectricalData(data);
+        busNumber++;
+    }
+
+    // Load
+    for(auto itlo = m_loadList.begin(); itlo != m_loadList.end(); itlo++) {
+        Load* load = *itlo;
+        if(load->IsOnline()) {
+            int n = ((Bus*)load->GetParentList()[0])->GetEletricalData().number;
+            LoadElectricalData data = load->GetPUElectricalData(systemPowerBase);
+            if(data.loadType == CONST_IMPEDANCE)
+                yBus[n][n] += std::complex<double>(data.activePower, -data.reactivePower);
+        }
+    }
+
+    // Capacitor
+    for(auto itca = m_capacitorList.begin(); itca != m_capacitorList.end(); itca++) {
+        Capacitor* capacitor = *itca;
+        if(capacitor->IsOnline()) {
+            int n = ((Bus*)capacitor->GetParentList()[0])->GetEletricalData().number;
+            CapacitorElectricalData data = capacitor->GetPUElectricalData(systemPowerBase);
+            yBus[n][n] += std::complex<double>(0.0, data.reactivePower);
+        }
+    }
+
+    // Inductor
+    for(auto itin = m_inductorList.begin(); itin != m_inductorList.end(); itin++) {
+        Inductor* inductor = *itin;
+        if(inductor->IsOnline()) {
+            int n = ((Bus*)inductor->GetParentList()[0])->GetEletricalData().number;
+            InductorElectricalData data = inductor->GetPUElectricalData(systemPowerBase);
+            yBus[n][n] += std::complex<double>(0.0, -data.reactivePower);
+        }
+    }
+
+    // Power line
+    for(auto itl = m_lineList.begin(); itl != m_lineList.end(); itl++) {
+        Line* line = *itl;
+        if(line->IsOnline()) {
+            LineElectricalData data = line->GetElectricalData();
+
+            int n1 = ((Bus*)line->GetParentList()[0])->GetEletricalData().number;
+            int n2 = ((Bus*)line->GetParentList()[1])->GetEletricalData().number;
+
+            yBus[n1][n2] -= 1.0 / std::complex<double>(data.resistance, data.indReactance);
+            yBus[n2][n1] -= 1.0 / std::complex<double>(data.resistance, data.indReactance);
+
+            yBus[n1][n1] += 1.0 / std::complex<double>(data.resistance, data.indReactance);
+            yBus[n2][n2] += 1.0 / std::complex<double>(data.resistance, data.indReactance);
+
+            yBus[n1][n1] += std::complex<double>(0.0, data.capSusceptance / 2.0);
+            yBus[n2][n2] += std::complex<double>(0.0, data.capSusceptance / 2.0);
+        }
+    }
+
+    // Transformer
+    for(auto itt = m_transformerList.begin(); itt != m_transformerList.end(); ++itt) {
+        Transformer* transformer = *itt;
+
+        if(transformer->IsOnline()) {
+            TransformerElectricalData data = transformer->GetElectricalData();
+
+            int n1 = ((Bus*)transformer->GetParentList()[0])->GetEletricalData().number;
+            int n2 = ((Bus*)transformer->GetParentList()[1])->GetEletricalData().number;
+
+            // If the transformer have nominal turns ratio (1.0) and no phase shifting, it will be modelled as series
+            // impedance.
+            if(data.turnsRatio == 1.0 && data.phaseShift == 0.0) {
+                yBus[n1][n2] += -1.0 / std::complex<double>(data.resistance, data.indReactance);
+                yBus[n2][n1] += -1.0 / std::complex<double>(data.resistance, data.indReactance);
+
+                yBus[n1][n1] += 1.0 / std::complex<double>(data.resistance, data.indReactance);
+                yBus[n2][n2] += 1.0 / std::complex<double>(data.resistance, data.indReactance);
+            }
+            // If the transformer have no-nominal turn ratio and/or phase shifting, it will be modelled in a different
+            // way (see references).
+            //[Ref. 1: Elementos de analise de sistemas de potencia - Stevenson - pg. 232]
+            //[Ref. 2: http://www.ee.washington.edu/research/real/Library/Reports/Tap_Adjustments_in_AC_Load_Flows.pdf]
+            //[Ref. 3: http://www.columbia.edu/~dano/courses/power/notes/power/andersson1.pdf]
+            else {
+                // Complex turns ratio
+                double radPhaseShift = wxDegToRad(data.phaseShift);
+                std::complex<double> a = std::complex<double>(data.turnsRatio * std::cos(radPhaseShift),
+                                                              -data.turnsRatio * std::sin(radPhaseShift));
+
+                // Transformer admitance
+                std::complex<double> y = 1.0 / std::complex<double>(data.resistance, data.indReactance);
+
+                yBus[n1][n1] += y / std::pow(std::abs(a), 2.0);
+                yBus[n1][n2] += -(y / std::conj(a));
+                yBus[n2][n1] += -(y / a);
+                yBus[n2][n2] += y;
+            }
+        }
+    }
+
+    return true;
+}
+
+void ElectricCalculation::UpdateElementsPowerFlow(std::vector<std::complex<double> > voltage,
+                                                  std::vector<std::complex<double> > power,
+                                                  std::vector<BusType> busType,
+                                                  std::vector<ReactiveLimits> reactiveLimit,
+                                                  double systemPowerBase)
+{
+    // Buses voltages
+    for(int i = 0; i < (int)m_busList.size(); i++) {
+        Bus* bus = m_busList[i];
+        BusElectricalData data = bus->GetEletricalData();
+        data.voltage = voltage[i];
+        bus->SetElectricalData(data);
+    }
+
+    // Power line
+    for(int i = 0; i < (int)m_lineList.size(); i++) {
+        Line* line = m_lineList[i];
+        if(line->IsOnline()) {
+            int n1 = ((Bus*)line->GetParentList()[0])->GetEletricalData().number;
+            int n2 = ((Bus*)line->GetParentList()[1])->GetEletricalData().number;
+
+            LineElectricalData data = line->GetElectricalData();
+            std::complex<double> v1 = voltage[n1];
+            std::complex<double> v2 = voltage[n2];
+
+            data.current[0] = (v1 - v2) / std::complex<double>(data.resistance, data.indReactance) +
+                              v1 * std::complex<double>(0.0, data.capSusceptance / 2.0);
+            data.current[1] = (v2 - v1) / std::complex<double>(data.resistance, data.indReactance) +
+                              v2 * std::complex<double>(0.0, data.capSusceptance / 2.0);
+
+            data.powerFlow[0] = v1 * std::conj(data.current[0]);
+            data.powerFlow[1] = v2 * std::conj(data.current[1]);
+
+            if(data.powerFlow[0].real() > data.powerFlow[1].real())
+                line->SetPowerFlowDirection(PF_BUS1_TO_BUS2);
+            else
+                line->SetPowerFlowDirection(PF_BUS2_TO_BUS1);
+
+            line->SetElectricalData(data);
+        }
+    }
+
+    // Transformer
+    for(int i = 0; i < (int)m_transformerList.size(); i++) {
+        Transformer* transformer = m_transformerList[i];
+        if(transformer->IsOnline()) {
+            TransformerElectricalData data = transformer->GetElectricalData();
+            int n1 = ((Bus*)transformer->GetParentList()[0])->GetEletricalData().number;
+            int n2 = ((Bus*)transformer->GetParentList()[1])->GetEletricalData().number;
+            std::complex<double> v1 = voltage[n1];  // Primary voltage
+            std::complex<double> v2 = voltage[n2];  // Secondary voltage
+
+            // Transformer admitance
+            std::complex<double> y = 1.0 / std::complex<double>(data.resistance, data.indReactance);
+
+            if(data.turnsRatio == 1.0 && data.phaseShift == 0.0) {
+                data.current[0] = (v1 - v2) * y;
+                data.current[1] = (v2 - v1) * y;
+            } else {
+                double radPS = wxDegToRad(data.phaseShift);
+                std::complex<double> a =
+                    std::complex<double>(data.turnsRatio * std::cos(radPS), -data.turnsRatio * std::sin(radPS));
+
+                data.current[0] = v1 * (y / std::pow(std::abs(a), 2)) - v2 * (y / std::conj(a));
+                data.current[1] = -v1 * (y / a) + v2 * y;
+            }
+
+            data.powerFlow[0] = v1 * std::conj(data.current[0]);
+            data.powerFlow[1] = v2 * std::conj(data.current[1]);
+
+            if(data.powerFlow[0].real() > data.powerFlow[1].real())
+                transformer->SetPowerFlowDirection(PF_BUS1_TO_BUS2);
+            else
+                transformer->SetPowerFlowDirection(PF_BUS2_TO_BUS1);
+
+            transformer->SetElectricaData(data);
+        }
+    }
+
+    // Synchronous machines
+    for(int i = 0; i < (int)m_busList.size(); i++) {
+        Bus* bus = m_busList[i];
+        BusElectricalData data = bus->GetEletricalData();
+
+        // Get the synchronous machines connected and calculate the load power on the bus.
+        std::vector<SyncGenerator*> syncGeneratorsOnBus;
+        std::vector<SyncMotor*> syncMotorsOnBus;
+        std::complex<double> loadPower(0.0, 0.0);
+
+        for(auto itsg = m_syncGeneratorList.begin(); itsg != m_syncGeneratorList.end(); itsg++) {
+            SyncGenerator* syncGenerator = *itsg;
+            if(bus == syncGenerator->GetParentList()[0] && syncGenerator->IsOnline())
+                syncGeneratorsOnBus.push_back(syncGenerator);
+        }
+        for(auto itsm = m_syncMotorList.begin(); itsm != m_syncMotorList.end(); itsm++) {
+            SyncMotor* syncMotor = *itsm;
+            if(bus == syncMotor->GetParentList()[0] && syncMotor->IsOnline()) {
+                syncMotorsOnBus.push_back(syncMotor);
+                SyncMotorElectricalData childData = syncMotor->GetPUElectricalData(systemPowerBase);
+                loadPower += std::complex<double>(childData.activePower, 0.0);
+            }
+        }
+        for(auto itlo = m_loadList.begin(); itlo != m_loadList.end(); itlo++) {
+            Load* load = *itlo;
+            if(bus == load->GetParentList()[0] && load->IsOnline()) {
+                LoadElectricalData childData = load->GetPUElectricalData(systemPowerBase);
+                if(childData.loadType == CONST_POWER)
+                    loadPower += std::complex<double>(childData.activePower, childData.reactivePower);
+
+                if(childData.activePower >= 0.0)
+                    load->SetPowerFlowDirection(PF_TO_ELEMENT);
+                else
+                    load->SetPowerFlowDirection(PF_TO_BUS);
+            }
+        }
+        for(auto itim = m_indMotorList.begin(); itim != m_indMotorList.end(); itim++) {
+            IndMotor* indMotor = *itim;
+            if(bus == indMotor->GetParentList()[0] && indMotor->IsOnline()) {
+                IndMotorElectricalData childData = indMotor->GetPUElectricalData(systemPowerBase);
+                loadPower += std::complex<double>(childData.activePower, childData.reactivePower);
+
+                if(childData.activePower >= 0.0)
+                    indMotor->SetPowerFlowDirection(PF_TO_ELEMENT);
+                else
+                    indMotor->SetPowerFlowDirection(PF_TO_BUS);
+            }
+        }
+
+        // Set the sync generator power
+        for(auto itsg = syncGeneratorsOnBus.begin(); itsg != syncGeneratorsOnBus.end(); itsg++) {
+            SyncGenerator* generator = *itsg;
+            if(generator->IsOnline()) {
+                SyncGeneratorElectricalData childData = generator->GetElectricalData();
+
+                if(busType[i] == BUS_SLACK) {
+                    double activePower =
+                        (power[i].real() + loadPower.real()) * systemPowerBase / (double)(syncGeneratorsOnBus.size());
+
+                    switch(childData.activePowerUnit) {
+                        case UNIT_PU: {
+                            activePower /= systemPowerBase;
+                        } break;
+                        case UNIT_kW: {
+                            activePower /= 1e3;
+                        } break;
+                        case UNIT_MW: {
+                            activePower /= 1e6;
+                        } break;
+                        default:
+                            break;
+                    }
+
+                    childData.activePower = activePower;
+                }
+                if(busType[i] == BUS_PV || busType[i] == BUS_SLACK) {
+                    // double reactivePower = (power[i].imag() + loadPower.imag()) * systemPowerBase /
+                    //                       (double)(syncGeneratorsOnBus.size() + syncMotorsOnBus.size());
+                    SyncGeneratorElectricalData childData_PU = generator->GetPUElectricalData(systemPowerBase);
+
+                    double reactivePower = (power[i].imag() + loadPower.imag()) * systemPowerBase;
+
+                    if(reactiveLimit[i].limitReached == RL_MAX_REACHED)
+                        reactivePower *= (childData_PU.maxReactive / reactiveLimit[i].maxLimit);
+
+                    else if(reactiveLimit[i].limitReached == RL_MIN_REACHED)
+                        reactivePower *= (childData_PU.minReactive / reactiveLimit[i].minLimit);
+
+                    else
+                        reactivePower /= (double)(syncGeneratorsOnBus.size() + syncMotorsOnBus.size());
+
+                    switch(childData.reactivePowerUnit) {
+                        case UNIT_PU: {
+                            reactivePower /= systemPowerBase;
+                        } break;
+                        case UNIT_kVAr: {
+                            reactivePower /= 1e3;
+                        } break;
+                        case UNIT_MVAr: {
+                            reactivePower /= 1e6;
+                        } break;
+                        default:
+                            break;
+                    }
+                    childData.reactivePower = reactivePower;
+                }
+
+                if(childData.activePower >= 0.0)
+                    generator->SetPowerFlowDirection(PF_TO_BUS);
+                else
+                    generator->SetPowerFlowDirection(PF_TO_ELEMENT);
+
+                generator->SetElectricalData(childData);
+            }
+        }
+
+        // Set the sync motor reactive power
+        double exceededReactive = 0.0;
+        int numMachines = syncGeneratorsOnBus.size() + syncMotorsOnBus.size();
+        for(auto itsm = syncMotorsOnBus.begin(); itsm != syncMotorsOnBus.end(); itsm++) {
+            SyncMotor* syncMotor = *itsm;
+            SyncMotorElectricalData childData = syncMotor->GetElectricalData();
+
+            bool reachedMachineLimit = false;
+
+            if(busType[i] == BUS_PV || busType[i] == BUS_SLACK) {
+                // double reactivePower = (power[i].imag() + loadPower.imag()) * systemPowerBase /
+                //                       (double)(syncGeneratorsOnBus.size() + syncMotorsOnBus.size());
+
+                SyncMotorElectricalData childData_PU = syncMotor->GetPUElectricalData(systemPowerBase);
+
+                double reactivePower = power[i].imag() + loadPower.imag();
+
+                // Bus reachd maximum reactive limit.
+                if(reactiveLimit[i].limitReached == RL_MAX_REACHED)
+                    reactivePower *= (childData_PU.maxReactive / reactiveLimit[i].maxLimit);
+                // Bus reached minimum reactive limit.
+                else if(reactiveLimit[i].limitReached == RL_MIN_REACHED)
+                    reactivePower *= (childData_PU.minReactive / reactiveLimit[i].minLimit);
+                // Bus didn't reach any limits
+                else {
+                    reactivePower /= (double)(numMachines);
+                    if(childData_PU.haveMaxReactive && (reactivePower > childData_PU.maxReactive)) {
+                        exceededReactive += reactivePower - childData_PU.maxReactive;
+                        reactivePower = childData_PU.maxReactive;
+                        reachedMachineLimit = true;
+                    } else if(childData_PU.haveMinReactive && (reactivePower < childData_PU.minReactive)) {
+                        exceededReactive += reactivePower - childData_PU.minReactive;
+                        reactivePower = childData_PU.minReactive;
+                        reachedMachineLimit = true;
+                    } else if((!childData_PU.haveMaxReactive && reactiveLimit[i].limitReached == RL_MAX_REACHED) ||
+                              (!childData_PU.haveMinReactive && reactiveLimit[i].limitReached == RL_MIN_REACHED) ||
+                              (!childData_PU.haveMaxReactive && !childData_PU.haveMaxReactive)) {
+                        reactivePower += exceededReactive;
+                        exceededReactive = 0.0;
+                    }
+                }
+
+                reactivePower *= systemPowerBase;
+
+                switch(childData.reactivePowerUnit) {
+                    case UNIT_PU: {
+                        reactivePower /= systemPowerBase;
+                    } break;
+                    case UNIT_kVAr: {
+                        reactivePower /= 1e3;
+                    } break;
+                    case UNIT_MVAr: {
+                        reactivePower /= 1e6;
+                    } break;
+                    default:
+                        break;
+                }
+                childData.reactivePower = reactivePower;
+            }
+
+            if(childData.activePower > 0.0)
+                syncMotor->SetPowerFlowDirection(PF_TO_ELEMENT);
+            else
+                syncMotor->SetPowerFlowDirection(PF_TO_BUS);
+
+            syncMotor->SetElectricalData(childData);
+
+            if(reachedMachineLimit) {
+                syncMotorsOnBus.erase(itsm);
+                itsm = syncMotorsOnBus.begin();
+            }
+        }
+    }
+}