Base class for electrical calculations providing general utility methods. More...

#include <ElectricCalculation.h>

Inheritance diagram for ElectricCalculation:

[legend]

Public Member Functions
	ElectricCalculation ()
	Constructor.

	~ElectricCalculation ()
	Destructor.

virtual void	GetElementsFromList (std::vector< Element * > elementList)
	Separate the power elements from a generic list.

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).

virtual bool	InvertMatrix (std::vector< std::vector< std::complex< double > > > matrix, std::vector< std::vector< std::complex< double > > > &inverse)
	Invert a matrix.

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.

void	ABCtoDQ0 (std::complex< double > complexValue, double angle, double &dValue, double &qValue)
	Convert a complex phasor in ABC representation to DQ components.

void	DQ0toABC (double dValue, double qValue, double angle, std::complex< double > &complexValue)
	Convert DQ components to a complex phasor in ABC representation.

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).

std::vector< double >	GaussianElimination (std::vector< std::vector< double > > matrix, std::vector< double > array)
	Solve a linear system using Gaussian elimination (real version).

Machines::SyncMachineModel	GetMachineModel (SyncGenerator *generator)
	Get the synchronous machine model used by the generator based on user-defined parameters.

std::vector< std::complex< double > >	ComplexMatrixTimesVector (std::vector< std::vector< std::complex< double > > > matrix, std::vector< std::complex< double > > vector)
	Multiply a complex matrix by a complex vector.

void	GetLUDecomposition (std::vector< std::vector< std::complex< double > > > matrix, std::vector< std::vector< std::complex< double > > > &matrixL, std::vector< std::vector< std::complex< double > > > &matrixU)
	Compute the LU decomposition of a matrix.

std::vector< std::complex< double > >	LUEvaluate (std::vector< std::vector< std::complex< double > > > u, std::vector< std::vector< std::complex< double > > > l, std::vector< std::complex< double > > b)
	Solve a linear system using LU decomposition.

bool	GetParentBus (Element childElement, Bus &parentBus)
	Get the parent bus of a given shunt element.

bool	GetParentBus (Element childElement, Bus &parentBus1, Bus *&parentBus2)
	Get the parent buses of a two-terminal element (branch).

bool	CalculateEMTElementsAdmittance (const double &basePower, wxString &errorMsg)
	Calculate the admittance of EMT elements.

bool	CalculateEMTElementsPower (const double &basePower, wxString &errorMsg, bool updateCurrent=true)
	Calculate the power of EMT elements.

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.

const std::vector< PowerElement * >	GetPowerElementList () const
	Get the power elements of the system (use GetElementsFromList first).

const std::vector< Bus * >	GetBusList () const
	Get the buses of the system (use GetElementsFromList first).

const std::vector< Capacitor * >	GetCapacitorList () const
	Get the capacitors of the system (use GetElementsFromList first).

const std::vector< IndMotor * >	GetIndMotorList () const
	Get the induction motors of the system (use GetElementsFromList first).

const std::vector< Inductor * >	GetInductorList () const
	Get the inductors of the system (use GetElementsFromList first).

const std::vector< Line * >	GetLineList () const
	Get the lines of the system (use GetElementsFromList first).

const std::vector< Load * >	GetLoadList () const
	Get the loads of the system (use GetElementsFromList first).

const std::vector< SyncGenerator * >	GetSyncGeneratorList () const
	Get the synchronous generators of the system (use GetElementsFromList first).

const std::vector< SyncMotor * >	GetSyncMotorList () const
	Get the synchronous motors of the system (use GetElementsFromList first).

const std::vector< Transformer * >	GetTransformerList () const
	Get the transformers of the system (use GetElementsFromList first).

const std::vector< HarmCurrent * >	GetHarmCurrentList () const
	Get the harmonic current source of the system (use GetElementsFromList first).

const std::vector< EMTElement * >	GetEMTElementList () const
	Get the electromagnetic element list of the system (use GetElementsFromList first).

Protected Member Functions
void	GetNextConnection (const unsigned int &checkBusNumber, const std::vector< std::vector< std::complex< double > > > &yBus, std::vector< bool > &connToSlack)
	Recursively check if a bus is electrically connected to the slack bus.

void	DistributeReactivePower (std::vector< ReactiveMachine > &machines, double qTotal)
	Distribute reactive power among synchronous machines connected to the same bus.

Protected Attributes
std::vector< PowerElement * >	m_powerElementList
	List of power elements in the system.

std::vector< Bus * >	m_busList
	List of buses in the system.

std::vector< Capacitor * >	m_capacitorList
	List of capacitor elements in the system.

std::vector< IndMotor * >	m_indMotorList
	List of induction motors in the system.

std::vector< Inductor * >	m_inductorList
	List of inductors in the system.

std::vector< Line * >	m_lineList
	List of transmission lines in the system.

std::vector< Load * >	m_loadList
	List of load elements in the system.

std::vector< SyncGenerator * >	m_syncGeneratorList
	List of synchronous generators in the system.

std::vector< SyncMotor * >	m_syncMotorList
	List of synchronous motors in the system.

std::vector< Transformer * >	m_transformerList
	List of transformers in the system.

std::vector< HarmCurrent * >	m_harmCurrentList
	List of harmonic current sources in the system.

std::vector< EMTElement * >	m_emtElementList
	List of electromagnetic transient (EMT) elements in the system.

Detailed Description

Base class for electrical calculations providing general utility methods.

Author: Thales Lima Oliveira

Date: 09/01/2017

This class provides common utilities used by different electrical calculations such as power flow, short-circuit analysis and EMT simulations.

Definition at line 115 of file ElectricCalculation.h.

Constructor & Destructor Documentation

◆ ElectricCalculation()

ElectricCalculation::ElectricCalculation ( )

Constructor.

Definition at line 25 of file ElectricCalculation.cpp.

25{}

◆ ~ElectricCalculation()

ElectricCalculation::~ElectricCalculation ( )

Destructor.

Definition at line 26 of file ElectricCalculation.cpp.

26{}

Member Function Documentation

◆ ABCtoDQ0()

void ElectricCalculation::ABCtoDQ0	(	std::complex< double >	complexValue,
		double	angle,
		double &	dValue,
		double &	qValue
	)

Convert a complex phasor in ABC representation to DQ components.

Parameters

complexValue	Complex phasor value in ABC reference frame.
angle	Electrical angle used for the transformation.
dValue	Resulting direct-axis component.
qValue	Resulting quadrature-axis component.

Definition at line 829 of file ElectricCalculation.cpp.

{
    dValue = -std::real(complexValue) * std::sin(angle) + std::imag(complexValue) * std::cos(angle);
    qValue = std::real(complexValue) * std::cos(angle) + std::imag(complexValue) * std::sin(angle);
}

◆ CalculateEMTElementsAdmittance()

bool ElectricCalculation::CalculateEMTElementsAdmittance	(	const double &	basePower,
		wxString &	errorMsg
	)

Calculate the admittance of EMT elements.

Parameters

basePower	System base power.
errorMsg	Error message in case of failure.

Returns: True if the calculation was successful.

Definition at line 1057 of file ElectricCalculation.cpp.

{
    for (EMTElement* emtElement : m_emtElementList) {
        if (!emtElement->IsOnline()) continue;
 
        emtElement->UpdateData();
        if (!emtElement->CalculateCurrent(errorMsg)) return false;
        auto data = emtElement->GetEMTElementData();
        double baseCurrent = basePower / (sqrt(3.0) * data.baseVoltage);
        std::complex<double> current = data.currHarmonics[1] / baseCurrent; // Fundamental frequency current
        std::complex<double> y0 = current / data.puVoltage;
        data.y0 = y0;
        data.power = std::complex<double>(0.0, 0.0);
        emtElement->SetEMTElementData(data);
    }
    return true;
}

◆ CalculateEMTElementsPower()

bool ElectricCalculation::CalculateEMTElementsPower	(	const double &	basePower,
		wxString &	errorMsg,
		bool	updateCurrent = `true`
	)

Calculate the power of EMT elements.

Parameters

basePower	System base power.
errorMsg	Error message in case of failure.
updateCurrent	If true, update the element currents during the calculation unsing EMTElement::CalculateCurrent(wxString&, const bool&).

Returns: True if the calculation was successful.

Definition at line 1075 of file ElectricCalculation.cpp.

{
    for (EMTElement* emtElement : m_emtElementList) {
        if (!emtElement->IsOnline()) continue;
 
        emtElement->UpdateData();
        if (updateCurrent) {
            if (!emtElement->CalculateCurrent(errorMsg)) return false;
        }
        auto data = emtElement->GetEMTElementData();
        double baseCurrent = basePower / (sqrt(3.0) * data.baseVoltage);
        //std::complex<double> current = data.currHarmonics[1];
        std::complex<double> current = data.currHarmonics[1] / baseCurrent; // Fundamental frequency current
        
        std::complex<double> i0 = data.y0 * data.puVoltage; // Current already injected via admittance matrix
        //std::complex<double> voltage = data.baseVoltage * data.puVoltage;
        //std::complex<double> power = (sqrt(3.0) * voltage * std::conj(current)) / basePower;
        std::complex<double> power = data.puVoltage * std::conj(current - i0);
        data.powerDiff = power - data.power;
        data.power = power;
        emtElement->SetEMTElementData(data);
 
        //wxString currentStr = data.name + ":\n";
        //for (auto const& [harm, current] : data.currHarmonics)
        //{
        //  currentStr += wxString::Format("I(%dh) = %f p.u.\n", harm, abs(current) / baseCurrent);
        //}
        //currentStr += wxString::Format("S = %e + j%e p.u.\n", data.power.real(), data.power.imag());
        //currentStr += wxString::Format("SDiff = %e p.u.\n", abs(data.powerDiff));
        //wxMessageBox(currentStr);
    }
    return true;
}

Here is the caller graph for this function:

◆ CalculateEMTPowerError()

double ElectricCalculation::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.

Parameters

voltage	Vector with bus voltages.
power	Vector with injected powers.
basePower	System base power.
errorMsg	Error message in case of failure.

Returns: The calculated power error.

Definition at line 1109 of file ElectricCalculation.cpp.

{
    // Update buses voltages
    for (Bus* bus : m_busList) {
        auto data = bus->GetElectricalData();
        if (data.isConnected) {
            data.voltage = voltage[data.number];
            bus->SetElectricalData(data);
        }
    }
 
    CalculateEMTElementsPower(basePower, errorMsg);
    // Update power array and calculate error
    double error = 0.0;
    for (EMTElement* emtElement : m_emtElementList) {
        if (!emtElement->IsOnline()) continue;
 
        if (!emtElement->GetParentList().empty()) {
            if (emtElement->GetParentList()[0] != nullptr) {
                int numBus = static_cast<Bus*>(emtElement->GetParentList()[0])->GetElectricalData().number;
                auto data = emtElement->GetEMTElementData();
                power[numBus] += data.powerDiff;
                if (error < std::abs(data.powerDiff))
                    error = std::abs(data.powerDiff);
            }
        }
    }
 
    return error;
}

Here is the call graph for this function:

◆ ComplexMatrixTimesVector()

std::vector< std::complex< double > > ElectricCalculation::ComplexMatrixTimesVector	(	std::vector< std::vector< std::complex< double > > >	matrix,
		std::vector< std::complex< double > >	vector
	)

Multiply a complex matrix by a complex vector.

Parameters

matrix	Complex matrix.
vector	Complex vector.

Returns: Resulting complex vector.

Definition at line 946 of file ElectricCalculation.cpp.

{
    std::vector<std::complex<double> > solution;
    for (unsigned int i = 0; i < matrix.size(); i++) {
        solution.push_back(std::complex<double>(0.0, 0.0));
 
        for (unsigned int j = 0; j < matrix.size(); j++) { solution[i] += matrix[i][j] * vector[j]; }
    }
 
    return solution;
}

◆ DistributeReactivePower()

void ElectricCalculation::DistributeReactivePower	(	std::vector< ReactiveMachine > &	machines,
		double	qTotal
	)

protected

Distribute reactive power among synchronous machines connected to the same bus.

The method ensures that reactive power limits of individual machines are respected during the distribution process.

Parameters

machines	List of machines participating in the reactive power sharing.
qTotal	Total reactive power to be distributed.

Definition at line 408 of file ElectricCalculation.cpp.

{
    int freeMachines = machines.size();
    double qRemaining = qTotal;
 
    bool changed = true;
 
    while (changed && freeMachines > 0) {
        changed = false;
 
        double qShare = qRemaining / freeMachines;
 
        for (auto& m : machines) {
            if (m.fixed)
                continue;
 
            if (m.hasMax && qShare > m.qMax) {
                m.q = m.qMax;
                m.fixed = true;
                qRemaining -= m.qMax;
                freeMachines--;
                changed = true;
            }
            else if (m.hasMin && qShare < m.qMin) {
                m.q = m.qMin;
                m.fixed = true;
                qRemaining -= m.qMin;
                freeMachines--;
                changed = true;
            }
            else {
                m.q = qShare;
            }
        }
    }
 
    if (freeMachines > 0) {
        double qShare = qRemaining / freeMachines;
        for (auto& m : machines) {
            if (!m.fixed)
                m.q = qShare;
        }
    }
}

Here is the caller graph for this function:

◆ DQ0toABC()

void ElectricCalculation::DQ0toABC	(	double	dValue,
		double	qValue,
		double	angle,
		std::complex< double > &	complexValue
	)

Convert DQ components to a complex phasor in ABC representation.

Parameters

dValue	Direct-axis component.
qValue	Quadrature-axis component.
angle	Electrical angle used for the transformation.
complexValue	Resulting complex phasor in ABC reference frame.

Definition at line 835 of file ElectricCalculation.cpp.

{
    double real = qValue * std::cos(angle) - dValue * std::sin(angle);
    double imag = qValue * std::sin(angle) + dValue * std::cos(angle);
    complexValue = std::complex<double>(real, imag);
}

◆ GaussianElimination() [1/2]

std::vector< double > ElectricCalculation::GaussianElimination	(	std::vector< std::vector< double > >	matrix,
		std::vector< double >	array
	)

Solve a linear system using Gaussian elimination (real version).

Parameters

matrix	Coefficient matrix of the system.
array	Right-hand side vector.

Returns: Solution vector of the linear system.

Definition at line 886 of file ElectricCalculation.cpp.

{
    //[Ref] https://en.wikipedia.org/wiki/Gaussian_elimination
 
    std::vector<double> solution;
 
    std::vector<std::vector<double> > triangMatrix;
    triangMatrix.resize(matrix.size());
    for (unsigned int i = 0; i < matrix.size(); i++) { triangMatrix[i].resize(matrix.size()); }
 
    for (unsigned int i = 0; i < matrix.size(); i++) { solution.push_back(array[i]); }
 
    for (unsigned int i = 0; i < matrix.size(); i++) {
        for (unsigned int j = 0; j < matrix.size(); j++) { triangMatrix[i][j] = matrix[i][j]; }
    }
 
    for (unsigned int k = 0; k < matrix.size(); k++) {
        unsigned int k1 = k + 1;
        for (unsigned int i = k; i < matrix.size(); i++) {
            if (triangMatrix[i][k] != 0.0) {
                for (unsigned int j = k1; j < matrix.size(); j++) {
                    triangMatrix[i][j] = triangMatrix[i][j] / triangMatrix[i][k];
                }
                solution[i] = solution[i] / triangMatrix[i][k];
            }
        }
        for (unsigned int i = k1; i < matrix.size(); i++) {
            if (triangMatrix[i][k] != 0.0) {
                for (unsigned int j = k1; j < matrix.size(); j++) { triangMatrix[i][j] -= triangMatrix[k][j]; }
                solution[i] -= solution[k];
            }
        }
    }
    for (int i = static_cast<int>(matrix.size()) - 2; i >= 0; i--) {
        for (int j = static_cast<int>(matrix.size()) - 1; j >= i + 1; j--) {
            solution[i] -= triangMatrix[i][j] * solution[j];
        }
    }
 
    return solution;
}

◆ GaussianElimination() [2/2]

std::vector< std::complex< double > > ElectricCalculation::GaussianElimination	(	std::vector< std::vector< std::complex< double > > >	matrix,
		std::vector< std::complex< double > >	array
	)

Solve a linear system using Gaussian elimination (complex version).

Parameters

matrix	Coefficient matrix of the system.
array	Right-hand side vector.

Returns: Solution vector of the linear system.

Definition at line 842 of file ElectricCalculation.cpp.

{
    //[Ref] https://en.wikipedia.org/wiki/Gaussian_elimination
 
    std::vector<std::complex<double> > solution;
 
    std::vector<std::vector<std::complex<double> > > triangMatrix;
    triangMatrix.resize(matrix.size());
    for (unsigned int i = 0; i < matrix.size(); i++) { triangMatrix[i].resize(matrix.size()); }
 
    for (unsigned int i = 0; i < matrix.size(); i++) { solution.push_back(array[i]); }
 
    for (unsigned int i = 0; i < matrix.size(); i++) {
        for (unsigned int j = 0; j < matrix.size(); j++) { triangMatrix[i][j] = matrix[i][j]; }
    }
 
    for (unsigned int k = 0; k < matrix.size(); k++) {
        unsigned int k1 = k + 1;
        for (unsigned int i = k; i < matrix.size(); i++) {
            if (triangMatrix[i][k] != std::complex<double>(0.0, 0.0)) {
                for (unsigned int j = k1; j < matrix.size(); j++) {
                    triangMatrix[i][j] = triangMatrix[i][j] / triangMatrix[i][k];
                }
                solution[i] = solution[i] / triangMatrix[i][k];
            }
        }
        for (unsigned int i = k1; i < matrix.size(); i++) {
            if (triangMatrix[i][k] != std::complex<double>(0.0, 0.0)) {
                for (unsigned int j = k1; j < matrix.size(); j++) { triangMatrix[i][j] -= triangMatrix[k][j]; }
                solution[i] -= solution[k];
            }
        }
    }
    for (int i = static_cast<int>(matrix.size()) - 2; i >= 0; i--) {
        for (int j = static_cast<int>(matrix.size()) - 1; j >= i + 1; j--) {
            solution[i] -= triangMatrix[i][j] * solution[j];
        }
    }
 
    return solution;
}

◆ GetBusList()

const std::vector< Bus * > ElectricCalculation::GetBusList ( ) const

inline

Get the buses of the system (use GetElementsFromList first).

Returns: A list of bus elements.

Definition at line 298 of file ElectricCalculation.h.

298{ return m_busList; }

◆ GetCapacitorList()

const std::vector< Capacitor * > ElectricCalculation::GetCapacitorList ( ) const

inline

Get the capacitors of the system (use GetElementsFromList first).

Returns: A list of capacitor elements.

Definition at line 303 of file ElectricCalculation.h.

303{ return m_capacitorList; }

ElectricCalculation::m_capacitorList

std::vector< Capacitor * > m_capacitorList

List of capacitor elements in the system.

Definition ElectricCalculation.h:372

◆ GetElementsFromList()

void ElectricCalculation::GetElementsFromList ( std::vector< Element * > elementList )

virtual

Separate the power elements from a generic list.

Parameters

elementList List of generic elements.

Definition at line 27 of file ElectricCalculation.cpp.

{
    m_powerElementList.clear();
    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();
    m_harmCurrentList.clear();
    m_emtElementList.clear();
    // TODO: Bad design?
    for (auto it = elementList.begin(); it != elementList.end(); it++) {
        Element* element = *it;
        m_powerElementList.push_back(static_cast<PowerElement*>(element));
 
        if (Bus* bus = dynamic_cast<Bus*>(element))
            m_busList.push_back(bus);
        else if (Capacitor* capacitor = dynamic_cast<Capacitor*>(element))
            m_capacitorList.push_back(capacitor);
        else if (IndMotor* indMotor = dynamic_cast<IndMotor*>(element))
            m_indMotorList.push_back(indMotor);
        else if (Inductor* inductor = dynamic_cast<Inductor*>(element))
            m_inductorList.push_back(inductor);
        else if (Line* line = dynamic_cast<Line*>(element))
            m_lineList.push_back(line);
        else if (Load* load = dynamic_cast<Load*>(element))
            m_loadList.push_back(load);
        else if (SyncGenerator* syncGenerator = dynamic_cast<SyncGenerator*>(element))
            m_syncGeneratorList.push_back(syncGenerator);
        else if (SyncMotor* syncMotor = dynamic_cast<SyncMotor*>(element))
            m_syncMotorList.push_back(syncMotor);
        else if (Transformer* transformer = dynamic_cast<Transformer*>(element))
            m_transformerList.push_back(transformer);
        else if (HarmCurrent* harmCurrent = dynamic_cast<HarmCurrent*>(element))
            m_harmCurrentList.push_back(harmCurrent);
        else if (EMTElement* emtElement = dynamic_cast<EMTElement*>(element))
            m_emtElementList.push_back(emtElement);
    }
 
    // Set buses numbers
    int busNumber = 0;
    for (auto itb = m_busList.begin(); itb != m_busList.end(); itb++) {
        Bus* bus = *itb;
        BusElectricalData data = bus->GetElectricalData();
        data.number = busNumber;
        bus->SetElectricalData(data);
        busNumber++;
    }
}

Here is the caller graph for this function:

◆ GetEMTElementList()

const std::vector< EMTElement * > ElectricCalculation::GetEMTElementList ( ) const

inline

Get the electromagnetic element list of the system (use GetElementsFromList first).

Returns: A list of electromagnetic elements.

Definition at line 348 of file ElectricCalculation.h.

348{ return m_emtElementList; }

◆ GetHarmCurrentList()

const std::vector< HarmCurrent * > ElectricCalculation::GetHarmCurrentList ( ) const

inline

Get the harmonic current source of the system (use GetElementsFromList first).

Returns: A list of harmonic current sources elements.

Definition at line 343 of file ElectricCalculation.h.

343{ return m_harmCurrentList; }

◆ GetIndMotorList()

const std::vector< IndMotor * > ElectricCalculation::GetIndMotorList ( ) const

inline

Get the induction motors of the system (use GetElementsFromList first).

Returns: A list of induction motor elements.

Definition at line 308 of file ElectricCalculation.h.

308{ return m_indMotorList; }

◆ GetInductorList()

const std::vector< Inductor * > ElectricCalculation::GetInductorList ( ) const

inline

Get the inductors of the system (use GetElementsFromList first).

Returns: A list of inductor elements.

Definition at line 313 of file ElectricCalculation.h.

313{ return m_inductorList; }

◆ GetLineList()

const std::vector< Line * > ElectricCalculation::GetLineList ( ) const

inline

Get the lines of the system (use GetElementsFromList first).

Returns: A list of line elements.

Definition at line 318 of file ElectricCalculation.h.

318{ return m_lineList; }

◆ GetLoadList()

const std::vector< Load * > ElectricCalculation::GetLoadList ( ) const

inline

Get the loads of the system (use GetElementsFromList first).

Returns: A list of load elements.

Definition at line 323 of file ElectricCalculation.h.

323{ return m_loadList; }

◆ GetLUDecomposition()

void ElectricCalculation::GetLUDecomposition	(	std::vector< std::vector< std::complex< double > > >	matrix,
		std::vector< std::vector< std::complex< double > > > &	matrixL,
		std::vector< std::vector< std::complex< double > > > &	matrixU
	)

Compute the LU decomposition of a matrix.

Parameters

matrix	Matrix to be decomposed.
matrixL	Lower triangular matrix.
matrixU	Upper triangular matrix.

Definition at line 960 of file ElectricCalculation.cpp.

{
    // Doolittle method
    // [Ref] http://www3.nd.edu/~zxu2/acms40390F11/Alg-LU-Crout.pdf
    // [Ref] http://www.engr.colostate.edu/~thompson/hPage/CourseMat/Tutorials/CompMethods/doolittle.pdf
 
    int size = static_cast<int>(matrix.size());  // Decomposed matrix size.
 
    // Set upper and lower matrices sizes.
    matrixL.resize(size);
    matrixU.resize(size);
    for (int i = 0; i < size; i++) {
        matrixL[i].resize(size);
        matrixU[i].resize(size);
    }
 
    // First row of upper matrix and first column of lower matrix.
    for (int i = 0; i < size; i++) {
        matrixU[0][i] = matrix[0][i];
        matrixL[i][0] = matrix[i][0] / matrixU[0][0];
    }
 
    // Lower matrix main diagonal.
    for (int i = 1; i < size; i++) { matrixL[i][i] = std::complex<double>(1.0, 0.0); }
 
    for (int i = 1; i < size - 1; i++) {
        // Upper matrix main diagonal.
        matrixU[i][i] = matrix[i][i];
        for (int k = 0; k < i; k++) { matrixU[i][i] -= matrixL[i][k] * matrixU[k][i]; }
 
        // Others elements of upper matrix
        for (int j = i + 1; j < size; j++) {
            matrixU[i][j] = matrix[i][j];
            for (int k = 0; k < i; k++) { matrixU[i][j] -= matrixL[i][k] * matrixU[k][j]; }
        }
 
        // Lower matrix elements
        for (int j = i + 1; j < size; j++) {
            matrixL[j][i] = matrix[j][i];
            for (int k = 0; k < i; k++) { matrixL[j][i] -= matrixL[j][k] * matrixU[k][i]; }
            matrixL[j][i] = matrixL[j][i] / matrixU[i][i];
        }
    }
 
    // Last element of upper matrix.
    matrixU[size - 1][size - 1] = matrix[size - 1][size - 1];
    for (int k = 0; k < size - 1; k++) { matrixU[size - 1][size - 1] -= matrixL[size - 1][k] * matrixU[k][size - 1]; }
}

◆ GetMachineModel()

Machines::SyncMachineModel ElectricCalculation::GetMachineModel ( SyncGenerator * generator )

Get the synchronous machine model used by the generator based on user-defined parameters.

Parameters

generator Pointer to the synchronous generator.

Returns: Enumeration indicating the synchronous machine model.

Definition at line 929 of file ElectricCalculation.cpp.

{
    auto data = generator->GetElectricalData();
    if (data.transTd0 != 0.0) {
        if (data.transTq0 != 0.0) {
            if (data.subTd0 != 0.0 || data.subTq0 != 0.0) { return Machines::SM_MODEL_5; }
            return Machines::SM_MODEL_3;
        }
        else {
            if (data.subTd0 != 0.0 || data.subTq0 != 0.0) { return Machines::SM_MODEL_4; }
            return Machines::SM_MODEL_2;
        }
    }
 
    return Machines::SM_MODEL_1;
}

◆ GetNextConnection()

void ElectricCalculation::GetNextConnection	(	const unsigned int &	checkBusNumber,
		const std::vector< std::vector< std::complex< double > > > &	yBus,
		std::vector< bool > &	connToSlack
	)

protected

Recursively check if a bus is electrically connected to the slack bus.

Parameters

checkBusNumber	Index of the bus being checked.
yBus	Admittance matrix of the system.
connToSlack	Boolean vector indicating if each bus is connected to the slack bus.

Definition at line 392 of file ElectricCalculation.cpp.

{
    for (unsigned int i = 0; i < yBus.size(); i++) {
 
        if (i != checkBusNumber && !connToSlack[i]) {
            if (std::abs(yBus[i][checkBusNumber]) > 1e-6)
            {
                connToSlack[i] = true;
                GetNextConnection(i, yBus, connToSlack);
            }
        }
    }
}

Here is the call graph for this function:

Here is the caller graph for this function:

◆ GetParentBus() [1/2]

bool ElectricCalculation::GetParentBus	(	Element *	childElement,
		Bus *&	parentBus
	)

Get the parent bus of a given shunt element.

Parameters

childElement	Pointer to the element.
parentBus	Pointer to the parent bus.

Returns: True if the parent bus was successfully found.

Definition at line 1036 of file ElectricCalculation.cpp.

{
    auto parentList = childElement->GetParentList();
    if (parentList.size() < 1) return false;
    parentBus = static_cast<Bus*>(parentList[0]);
    if (parentBus == nullptr) return false;
    if (parentBus->GetElectricalData().number < 0) return false;
    return true;
}

Here is the call graph for this function:

◆ GetParentBus() [2/2]

bool ElectricCalculation::GetParentBus	(	Element *	childElement,
		Bus *&	parentBus1,
		Bus *&	parentBus2
	)

Get the parent buses of a two-terminal element (branch).

Parameters

childElement	Pointer to the element.
parentBus1	Pointer to the first parent bus.
parentBus2	Pointer to the second parent bus.

Returns: True if both parent buses were successfully found.

Definition at line 1046 of file ElectricCalculation.cpp.

{
    auto parentList = childElement->GetParentList();
    if (parentList.size() < 2) return false;
    parentBus1 = static_cast<Bus*>(parentList[0]);
    parentBus2 = static_cast<Bus*>(parentList[1]);
    if (parentBus1 == nullptr || parentBus2 == nullptr) return false;
    if (parentBus1->GetElectricalData().number < 0 || parentBus2->GetElectricalData().number < 0) return false;
    return true;
}

Here is the call graph for this function:

◆ GetPowerElementList()

const std::vector< PowerElement * > ElectricCalculation::GetPowerElementList ( ) const

inline

Get the power elements of the system (use GetElementsFromList first).

Returns: A list of power elements.

Definition at line 293 of file ElectricCalculation.h.

293{ return m_powerElementList; }

◆ GetSyncGeneratorList()

const std::vector< SyncGenerator * > ElectricCalculation::GetSyncGeneratorList ( ) const

inline

Get the synchronous generators of the system (use GetElementsFromList first).

Returns: A list of synchronous generator elements.

Definition at line 328 of file ElectricCalculation.h.

328{ return m_syncGeneratorList; }

◆ GetSyncMotorList()

const std::vector< SyncMotor * > ElectricCalculation::GetSyncMotorList ( ) const

inline

Get the synchronous motors of the system (use GetElementsFromList first).

Returns: A list of synchronous motor elements.

Definition at line 333 of file ElectricCalculation.h.

333{ return m_syncMotorList; }

◆ GetTransformerList()

const std::vector< Transformer * > ElectricCalculation::GetTransformerList ( ) const

inline

Get the transformers of the system (use GetElementsFromList first).

Returns: A list of transformer elements.

Definition at line 338 of file ElectricCalculation.h.

338{ return m_transformerList; }

◆ GetYBus()

bool ElectricCalculation::GetYBus	(	std::vector< std::vector< std::complex< double > > > &	yBus,
		double	systemPowerBase,
		YBusSequence	sequence = `POSITIVE_SEQ`,
		bool	includeSyncMachines = `false`,
		bool	allLoadsAsImpedances = `false`,
		bool	usePowerFlowVoltagesOnImpedances = `false`
	)

virtual

Get the admittance matrix from the list of elements (use GetElementsFromList first).

Parameters

yBus	Admittance matrix. The previous content will be erased.
systemPowerBase	Base power of the system.
sequence	Sequence of admittance matrix (positive, negative and zero).
includeSyncMachines	Include the synchronous machines on calculation.

Returns: Return true if was possible to build the admittance matrix.

Definition at line 81 of file ElectricCalculation.cpp.

{
    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 and reset connection status
    int busNumber = 0;
    for (auto itb = m_busList.begin(); itb != m_busList.end(); itb++) {
        Bus* bus = *itb;
        BusElectricalData data = bus->GetElectricalData();
        data.number = busNumber;
        data.isConnected = true;
        //data.voltage = std::complex<double>(1.0, 0.0);
        bus->SetElectricalData(data);
        busNumber++;
    }
 
    // Load
    for (auto it = m_loadList.begin(), itEnd = m_loadList.end(); it != itEnd; ++it) {
        Load* load = *it;
        if (load->IsOnline()) {
            int n = static_cast<Bus*>(load->GetParentList()[0])->GetElectricalData().number;
            LoadElectricalData data = load->GetPUElectricalData(systemPowerBase);
            if (data.loadType == CONST_IMPEDANCE || allLoadsAsImpedances) {
                std::complex<double> yLoad = std::complex<double>(data.activePower, -data.reactivePower);
                if (allLoadsAsImpedances && usePowerFlowVoltagesOnImpedances) {
                    std::complex<double> v = static_cast<Bus*>(load->GetParentList()[0])->GetElectricalData().voltage;
                    yLoad /= (std::abs(v) * std::abs(v));
                }
                yBus[n][n] += yLoad;
            }
        }
    }
 
    // Capacitor
    for (auto it = m_capacitorList.begin(), itEnd = m_capacitorList.end(); it != itEnd; ++it) {
        Capacitor* capacitor = *it;
        if (capacitor->IsOnline()) {
            int n = static_cast<Bus*>(capacitor->GetParentList()[0])->GetElectricalData().number;
            CapacitorElectricalData data = capacitor->GetPUElectricalData(systemPowerBase);
            yBus[n][n] += std::complex<double>(0.0, data.reactivePower);
        }
    }
 
    // Inductor
    for (auto it = m_inductorList.begin(), itEnd = m_inductorList.end(); it != itEnd; ++it) {
        Inductor* inductor = *it;
        if (inductor->IsOnline()) {
            int n = static_cast<Bus*>(inductor->GetParentList()[0])->GetElectricalData().number;
            InductorElectricalData data = inductor->GetPUElectricalData(systemPowerBase);
            yBus[n][n] += std::complex<double>(0.0, -data.reactivePower);
        }
    }
 
    //EMT Elements
    for (EMTElement* emtElement : m_emtElementList) {
        if (emtElement->IsOnline()) {
            auto data = emtElement->GetEMTElementData();
            int n = static_cast<Bus*>(emtElement->GetParentList()[0])->GetElectricalData().number;
            yBus[n][n] += data.y0;
            //wxMessageBox(wxString::Format("%f +j %f", data.y0.real(), data.y0.imag()));
        }
    }
 
    // Power line
    for (auto it = m_lineList.begin(), itEnd = m_lineList.end(); it != itEnd; ++it) {
        Line* line = *it;
        if (line->IsOnline()) {
            LineElectricalData data = line->GetPUElectricalData(systemPowerBase);
 
            int n1 = static_cast<Bus*>(line->GetParentList()[0])->GetElectricalData().number;
            int n2 = static_cast<Bus*>(line->GetParentList()[1])->GetElectricalData().number;
 
            switch (sequence) {
            case POSITIVE_SEQ:
            case NEGATIVE_SEQ: {
                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);
            } break;
            case ZERO_SEQ: {
                yBus[n1][n2] -= 1.0 / std::complex<double>(data.zeroResistance, data.zeroIndReactance);
                yBus[n2][n1] -= 1.0 / std::complex<double>(data.zeroResistance, data.zeroIndReactance);
 
                yBus[n1][n1] += 1.0 / std::complex<double>(data.zeroResistance, data.zeroIndReactance);
                yBus[n2][n2] += 1.0 / std::complex<double>(data.zeroResistance, data.zeroIndReactance);
 
                yBus[n1][n1] += std::complex<double>(0.0, data.zeroCapSusceptance / 2.0);
                yBus[n2][n2] += std::complex<double>(0.0, data.zeroCapSusceptance / 2.0);
            }
            }
        }
    }
 
    // Transformer
    for (auto it = m_transformerList.begin(), itEnd = m_transformerList.end(); it != itEnd; ++it) {
        Transformer* transformer = *it;
 
        if (transformer->IsOnline()) {
            TransformerElectricalData data = transformer->GetPUElectricalData(systemPowerBase);
 
            int n1 = static_cast<Bus*>(transformer->GetParentList()[0])->GetElectricalData().number;
            int n2 = static_cast<Bus*>(transformer->GetParentList()[1])->GetElectricalData().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 && sequence != ZERO_SEQ) {
                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 if (sequence != ZERO_SEQ) {
                // 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);
 
                if (sequence == POSITIVE_SEQ) {
                    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;
                }
                else if (sequence == NEGATIVE_SEQ) {
                    yBus[n1][n1] += y / std::pow(std::abs(a), 2.0);
                    yBus[n1][n2] += -(y / a);
                    yBus[n2][n1] += -(y / std::conj(a));
                    yBus[n2][n2] += y;
                }
            }
            else if (sequence == ZERO_SEQ) {
                switch (data.connection) {
                case GWYE_GWYE: {
                    std::complex<double> y =
                        1.0 /
                        std::complex<double>(
                            data.zeroResistance + 3.0 * (data.primaryGrndResistance + data.secondaryGrndResistance),
                            data.zeroIndReactance +
                            3.0 * (data.primaryGrndReactance + data.secondaryGrndReactance));
                    std::complex<double> a = std::complex<double>(data.turnsRatio, 0.0);
 
                    yBus[n1][n1] += y / (a * a);
                    yBus[n1][n2] += -(y / a);
                    yBus[n2][n1] += -(y / a);
                    yBus[n2][n2] += y;
                } break;
                case DELTA_GWYE: {
                    std::complex<double> y =
                        1.0 / std::complex<double>(data.zeroResistance + 3.0 * (data.secondaryGrndResistance),
                            data.zeroIndReactance + 3.0 * (data.secondaryGrndReactance));
                    yBus[n2][n2] += y;
                    yBus[n1][n1] += 1e-4;
                    yBus[n1][n2] += 1e-4;
                    yBus[n2][n1] += 1e-4;
                    break;
                }
                case GWYE_DELTA: {
                    std::complex<double> y =
                        1.0 / std::complex<double>(data.zeroResistance + 3.0 * (data.primaryGrndResistance),
                            data.zeroIndReactance + 3.0 * (data.primaryGrndReactance));
                    yBus[n1][n1] += y;
                    yBus[n2][n2] += 1e-4;
                    yBus[n1][n2] += 1e-4;
                    yBus[n2][n1] += 1e-4;
                    break;
                }
                default:
                    break;
                }
            }
        }
    }
 
    if (includeSyncMachines) {
        // Synchronous generator
        for (auto it = m_syncGeneratorList.begin(), itEnd = m_syncGeneratorList.end(); it != itEnd; ++it) {
            SyncGenerator* syncGenerator = *it;
            if (syncGenerator->IsOnline()) {
                int n = static_cast<Bus*>(syncGenerator->GetParentList()[0])->GetElectricalData().number;
                SyncGeneratorElectricalData data = syncGenerator->GetPUElectricalData(systemPowerBase);
                switch (sequence) {
                case POSITIVE_SEQ: {
                    yBus[n][n] += 1.0 / std::complex<double>(data.positiveResistance, data.positiveReactance);
                } break;
                case NEGATIVE_SEQ: {
                    yBus[n][n] += 1.0 / std::complex<double>(data.negativeResistance, data.negativeReactance);
                } break;
                case ZERO_SEQ: {
                    if (data.groundNeutral) {
                        yBus[n][n] += 1.0 / std::complex<double>(data.zeroResistance + 3.0 * data.groundResistance,
                            data.zeroReactance + 3.0 * data.groundReactance);
                    }
                } break;
                }
            }
        }
        // Synchronous motor
        for (auto it = m_syncMotorList.begin(), itEnd = m_syncMotorList.end(); it != itEnd; ++it) {
            SyncMotor* syncMotor = *it;
            if (syncMotor->IsOnline()) {
                int n = static_cast<Bus*>(syncMotor->GetParentList()[0])->GetElectricalData().number;
                SyncMotorElectricalData data = syncMotor->GetPUElectricalData(systemPowerBase);
                switch (sequence) {
                case POSITIVE_SEQ: {
                    yBus[n][n] += 1.0 / std::complex<double>(data.positiveResistance, data.positiveReactance);
                } break;
                case NEGATIVE_SEQ: {
                    yBus[n][n] += 1.0 / std::complex<double>(data.negativeResistance, data.negativeReactance);
                } break;
                case ZERO_SEQ: {
                    if (data.groundNeutral) {
                        yBus[n][n] += 1.0 / std::complex<double>(data.zeroResistance + 3.0 * data.groundResistance,
                            data.zeroReactance + 3.0 * data.groundReactance);
                    }
                } break;
                }
            }
        }
    }
 
    // Identify buses to slack bus (power flow)
    //if (sequence != ZERO_SEQ) {
    std::vector<bool> connectedToSlack(m_busList.size(), false);
    busNumber = 0;
    int slackBus = 0;
    std::vector<int> checkBusVector;
    for (auto* bus : m_busList)
    {
        if (bus->GetElectricalData().slackBus)
        {
            connectedToSlack[busNumber] = true;
            slackBus = busNumber;
            break;
        }
        busNumber++;
    }
    GetNextConnection(slackBus, yBus, connectedToSlack);
 
    bool hasIsolatedBus = false;
    for (auto isConnectedToSlack : connectedToSlack)
    {
        if (!isConnectedToSlack) hasIsolatedBus = true;
    }
 
    if (hasIsolatedBus) {
        // Update isolated and connected buses status and numbers
        busNumber = 0;
        for (unsigned int i = 0; i < m_busList.size(); ++i)
        {
            auto data = m_busList[i]->GetElectricalData();
            data.isConnected = connectedToSlack[i];
            if (connectedToSlack[i])
            {
                data.number = busNumber;
                busNumber++;
            }
            else {
                data.number = -1;
            }
            m_busList[i]->SetElectricalData(data);
        }
 
        // Create new Ybus without isolated buses (remove rows and columns)
        std::vector< std::vector< std::complex<double> > > newYBus;
        for (unsigned int i = 0; i < yBus.size(); ++i) {
            if (connectedToSlack[i]) {
                std::vector< std::complex<double> > newLine;
                for (unsigned int j = 0; j < yBus.size(); ++j) {
                    if (connectedToSlack[j]) {
                        newLine.push_back(yBus[i][j]);
                    }
                }
                newYBus.push_back(newLine);
            }
        }
 
        yBus.clear();
        yBus = newYBus;
    }
    //}
 
    return true;
}

Here is the call graph for this function:

Here is the caller graph for this function:

◆ InvertMatrix()

bool ElectricCalculation::InvertMatrix	(	std::vector< std::vector< std::complex< double > > >	matrix,
		std::vector< std::vector< std::complex< double > > > &	inverse
	)

virtual

Invert a matrix.

Parameters

matrix	Matrix to invert.
inverse	Inverted matrix. The previous content will be erased.

Returns: Return true if was possible to invert the matrix.

Definition at line 760 of file ElectricCalculation.cpp.

{
    int order = static_cast<int>(matrix.size());
 
    inverse.clear();
    // Fill the inverse matrix with identity.
    for (int i = 0; i < order; ++i) {
        std::vector<std::complex<double> > line;
        for (int j = 0; j < order; ++j) {
            line.push_back(i == j ? std::complex<double>(1.0, 0.0) : std::complex<double>(0.0, 0.0));
        }
        inverse.push_back(line);
    }
 
    // Check if a main diagonal value of the matrix is zero, if one is zero, try a linear combination to remove it.
    for (int i = 0; i < order; ++i) {
        for (int j = 0; j < order; ++j) {
            if (i == j && matrix[i][j] == std::complex<double>(0.0, 0.0)) {
                int row = 0;
                while (row < order) {
                    if (matrix[row][j] != std::complex<double>(0.0, 0.0)) {
                        for (int k = 0; k < order; ++k) {
                            matrix[i][k] += matrix[row][k];
                            inverse[i][k] += inverse[row][k];
                        }
                        break;
                    }
                    row++;
                }
                // If all line values are zero, the matrix is singular and the solution is impossible.
                if (row == order) return false;
            }
        }
    }
 
    // Linear combinations are made in both matrices, the goal is the input matrix become the identity. The final result
    // have two matrices: the identity and the inverse of the input.
    for (int i = 0; i < order; ++i) {
        for (int j = 0; j < order; ++j) {
            if (i != j) {
                if (matrix[i][i] == std::complex<double>(0.0, 0.0)) return false;
 
                std::complex<double> factor = matrix[j][i] / matrix[i][i];
                for (int k = 0; k < order; ++k) {
                    matrix[j][k] -= factor * matrix[i][k];
                    inverse[j][k] -= factor * inverse[i][k];
                }
            }
        }
    }
    // Main diagonal calculation.
    for (int i = 0; i < order; ++i) {
        for (int j = 0; j < order; ++j) {
            if (i == j) {
                if (matrix[i][j] == std::complex<double>(0.0, 0.0)) return false;
 
                std::complex<double> factor = (matrix[i][j] - std::complex<double>(1.0, 0.0)) / matrix[i][j];
                for (int k = 0; k < order; ++k) {
                    matrix[j][k] -= factor * matrix[i][k];
                    inverse[j][k] -= factor * inverse[i][k];
                }
            }
        }
    }
 
    return true;
}

Here is the caller graph for this function:

◆ LUEvaluate()

std::vector< std::complex< double > > ElectricCalculation::LUEvaluate	(	std::vector< std::vector< std::complex< double > > >	u,
		std::vector< std::vector< std::complex< double > > >	l,
		std::vector< std::complex< double > >	b
	)

Solve a linear system using LU decomposition.

Parameters

u	Upper triangular matrix.
l	Lower triangular matrix.
b	Right-hand side vector.

Returns: Solution vector.

Definition at line 1011 of file ElectricCalculation.cpp.

{
    int size = static_cast<int>(b.size());
    std::vector<std::complex<double> > x;
    std::vector<std::complex<double> > y;
    x.resize(size);
    y.resize(size);
 
    // Forward
    for (int i = 0; i < size; i++) {
        y[i] = b[i];
        for (int j = 0; j < i; j++) { y[i] -= l[i][j] * y[j]; }
        y[i] /= l[i][i];
    }
    // Backward
    for (int i = size - 1; i >= 0; i--) {
        x[i] = y[i];
        for (int j = i + 1; j < size; j++) { x[i] -= u[i][j] * x[j]; }
        x[i] /= u[i][i];
    }
    return x;
}

◆ UpdateElementsPowerFlow()

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
	)

virtual

Update the elements after the power flow calculation.

Parameters

voltage	Array with the buses voltages.
power	Array with the buses injected power.
busType	Array with the buses type.
reactiveLimit	Array with the reactive limit data.
systemPowerBase	Base power of the system.

Definition at line 453 of file ElectricCalculation.cpp.

{
    double zeroLimit = 1e-4;
    for (unsigned int i = 0; i < reactiveLimit.size(); ++i) {
        if (reactiveLimit[i].maxLimit > -zeroLimit && reactiveLimit[i].maxLimit < zeroLimit)
            reactiveLimit[i].maxLimit = zeroLimit;
        if (reactiveLimit[i].minLimit > -zeroLimit && reactiveLimit[i].minLimit < zeroLimit)
            reactiveLimit[i].minLimit = zeroLimit;
    }
    for (unsigned int i = 0; i < power.size(); ++i) {
        if (std::real(power[i]) > -zeroLimit && std::real(power[i]) < zeroLimit)
            power[i] = std::complex<double>(0.0, std::imag(power[i]));
        if (std::imag(power[i]) > -zeroLimit && std::imag(power[i]) < zeroLimit)
            power[i] = std::complex<double>(std::real(power[i]), 0.0);
    }
    // Buses
    for (int i = 0; i < (int)m_busList.size(); i++) {
        Bus* bus = m_busList[i];
        BusElectricalData data = bus->GetElectricalData();
        if (data.isConnected) {
            data.voltage = voltage[data.number];
            data.power = power[data.number];
            data.busType = busType[data.number];
        }
        else {
            data.voltage = std::complex<double>(0.0, 0.0);
            data.power = std::complex<double>(0.0, 0.0);
        }
        bus->SetElectricalData(data);
    }
 
    // Power line
    for (int i = 0; i < (int)m_lineList.size(); i++) {
        Line* line = m_lineList[i];
        if (line->IsOnline()) {
            auto dataBus1 = static_cast<Bus*>(line->GetParentList()[0])->GetElectricalData();
            auto dataBus2 = static_cast<Bus*>(line->GetParentList()[1])->GetElectricalData();
 
            if (dataBus1.isConnected && dataBus2.isConnected) {
                int n1 = dataBus1.number;
                int n2 = dataBus2.number;
 
                LineElectricalData data = line->GetElectricalData();
                LineElectricalData dataPU = line->GetPUElectricalData(systemPowerBase);
                std::complex<double> v1 = voltage[n1];
                std::complex<double> v2 = voltage[n2];
 
                data.current[0] = (v1 - v2) / std::complex<double>(dataPU.resistance, dataPU.indReactance) +
                    v1 * std::complex<double>(0.0, dataPU.capSusceptance / 2.0);
                data.current[1] = (v2 - v1) / std::complex<double>(dataPU.resistance, dataPU.indReactance) +
                    v2 * std::complex<double>(0.0, dataPU.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(PowerFlowDirection::PF_BUS1_TO_BUS2);
                else
                    line->SetPowerFlowDirection(PowerFlowDirection::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()) {
            auto dataBus1 = static_cast<Bus*>(transformer->GetParentList()[0])->GetElectricalData();
            auto dataBus2 = static_cast<Bus*>(transformer->GetParentList()[1])->GetElectricalData();
 
            if (dataBus1.isConnected && dataBus2.isConnected) {
                TransformerElectricalData data = transformer->GetElectricalData();
                TransformerElectricalData dataPU = transformer->GetPUElectricalData(systemPowerBase);
 
                int n1 = dataBus1.number;
                int n2 = dataBus2.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>(dataPU.resistance, dataPU.indReactance);
 
                if (dataPU.turnsRatio == 1.0 && dataPU.phaseShift == 0.0) {
                    data.current[0] = (v1 - v2) * y;
                    data.current[1] = (v2 - v1) * y;
                }
                else {
                    double radPS = wxDegToRad(dataPU.phaseShift);
                    std::complex<double> a =
                        std::complex<double>(dataPU.turnsRatio * std::cos(radPS), -dataPU.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(PowerFlowDirection::PF_BUS1_TO_BUS2);
                else
                    transformer->SetPowerFlowDirection(PowerFlowDirection::PF_BUS2_TO_BUS1);
 
                transformer->SetElectricaData(data);
            }
        }
    }
 
    // Ind motor
    for (unsigned int i = 0; i < m_indMotorList.size(); i++) {
        IndMotor* motor = m_indMotorList[i];
        auto data = motor->GetElectricalData();
        if (motor->IsOnline() && data.calcQInPowerFlow) {
            double reactivePower = data.qValue * systemPowerBase;
 
            switch (data.reactivePowerUnit) {
            case ElectricalUnit::UNIT_PU: {
                reactivePower /= systemPowerBase;
            } break;
            case ElectricalUnit::UNIT_kvar: {
                reactivePower /= 1e3;
            } break;
            case ElectricalUnit::UNIT_Mvar: {
                reactivePower /= 1e6;
            } break;
            default:
                break;
            }
 
            data.reactivePower = reactivePower;
 
            motor->SetElectricalData(data);
        }
    }
 
    //EMT Elements
    for (EMTElement* emtElement : m_emtElementList) {
 
        if (!emtElement->IsOnline()) {
            auto data = emtElement->GetEMTElementData();
            data.power = std::complex<double>(0.0, 0.0);
            data.currHarmonics[1] = std::complex<double>(0.0, 0.0);
            emtElement->SetEMTElementData(data);
        }
        else {
            emtElement->UpdateData();
            auto data = emtElement->GetEMTElementData();
            double baseCurrent = systemPowerBase / (sqrt(3.0) * data.baseVoltage);
            std::complex<double> current = data.currHarmonics[1] / baseCurrent; // Fundamental frequency current
            std::complex<double> power = data.puVoltage * std::conj(current);
            data.power = power;
            emtElement->SetEMTElementData(data);
        }
    }
 
    // Synchronous machines
    for (auto* bus : m_busList) {
        BusElectricalData data = bus->GetElectricalData();
        int i = data.number;
 
        if (data.isConnected) {
            // 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(PowerFlowDirection::PF_TO_ELEMENT);
                    else
                        load->SetPowerFlowDirection(PowerFlowDirection::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(PowerFlowDirection::PF_TO_ELEMENT);
                    else
                        indMotor->SetPowerFlowDirection(PowerFlowDirection::PF_TO_BUS);
                }
            }
 
            // Setup the reactive machines vector for the distribution of reactive power
            std::vector<ReactiveMachine> machines;
            for (auto* gen : syncGeneratorsOnBus) {
 
                auto data = gen->GetPUElectricalData(systemPowerBase);
 
                ReactiveMachine m;
                m.qMax = data.maxReactive;
                m.qMin = data.minReactive;
 
                m.hasMax = data.haveMaxReactive;
                m.hasMin = data.haveMinReactive;
 
                m.machine = gen;
                m.isGenerator = true;
 
                machines.push_back(m);
            }
            for (auto* mot : syncMotorsOnBus) {
 
                auto data = mot->GetPUElectricalData(systemPowerBase);
 
                ReactiveMachine m;
                m.qMax = data.maxReactive;
                m.qMin = data.minReactive;
 
                m.hasMax = data.haveMaxReactive;
                m.hasMin = data.haveMinReactive;
 
                m.machine = mot;
                m.isGenerator = false;
 
                machines.push_back(m);
            }
 
            // Distribute the reactive power among the synchronous machines and set their reactive power
            double qTotal = power[i].imag() + loadPower.imag();
            DistributeReactivePower(machines, qTotal);
 
            // Set the reactive power on the machines
            for (auto& m : machines) {
 
                if (m.isGenerator) {
 
                    SyncGenerator* gen = static_cast<SyncGenerator*>(m.machine);
                    auto data = gen->GetElectricalData();
 
                    double reactivePower = m.q * systemPowerBase;
 
                    switch (data.reactivePowerUnit) {
                    case ElectricalUnit::UNIT_PU:
                        reactivePower /= systemPowerBase;
                        break;
                    case ElectricalUnit::UNIT_kvar:
                        reactivePower /= 1e3;
                        break;
                    case ElectricalUnit::UNIT_Mvar:
                        reactivePower /= 1e6;
                        break;
                    default:
                        break;
                    }
 
                    data.reactivePower = reactivePower;
                    gen->SetElectricalData(data);
                }
                else {
 
                    SyncMotor* mot = static_cast<SyncMotor*>(m.machine);
                    auto data = mot->GetElectricalData();
 
                    double reactivePower = m.q * systemPowerBase;
 
                    switch (data.reactivePowerUnit) {
                    case ElectricalUnit::UNIT_PU:
                        reactivePower /= systemPowerBase;
                        break;
                    case ElectricalUnit::UNIT_kvar:
                        reactivePower /= 1e3;
                        break;
                    case ElectricalUnit::UNIT_Mvar:
                        reactivePower /= 1e6;
                        break;
                    default:
                        break;
                    }
 
                    data.reactivePower = reactivePower;
                    mot->SetElectricalData(data);
                }
            }
        }
    }
}

Here is the call graph for this function:

Member Data Documentation

◆ m_busList

std::vector<Bus*> ElectricCalculation::m_busList

protected

List of buses in the system.

Definition at line 370 of file ElectricCalculation.h.

◆ m_capacitorList

std::vector<Capacitor*> ElectricCalculation::m_capacitorList

protected

List of capacitor elements in the system.

Definition at line 372 of file ElectricCalculation.h.

◆ m_emtElementList

std::vector<EMTElement*> ElectricCalculation::m_emtElementList

protected

List of electromagnetic transient (EMT) elements in the system.

Definition at line 390 of file ElectricCalculation.h.

◆ m_harmCurrentList

std::vector<HarmCurrent*> ElectricCalculation::m_harmCurrentList

protected

List of harmonic current sources in the system.

Definition at line 388 of file ElectricCalculation.h.

◆ m_indMotorList

std::vector<IndMotor*> ElectricCalculation::m_indMotorList

protected

List of induction motors in the system.

Definition at line 374 of file ElectricCalculation.h.

◆ m_inductorList

std::vector<Inductor*> ElectricCalculation::m_inductorList

protected

List of inductors in the system.

Definition at line 376 of file ElectricCalculation.h.

◆ m_lineList

std::vector<Line*> ElectricCalculation::m_lineList

protected

List of transmission lines in the system.

Definition at line 378 of file ElectricCalculation.h.

◆ m_loadList

std::vector<Load*> ElectricCalculation::m_loadList

protected

List of load elements in the system.

Definition at line 380 of file ElectricCalculation.h.

◆ m_powerElementList

std::vector<PowerElement*> ElectricCalculation::m_powerElementList

protected

List of power elements in the system.

Definition at line 368 of file ElectricCalculation.h.

◆ m_syncGeneratorList

std::vector<SyncGenerator*> ElectricCalculation::m_syncGeneratorList

protected

List of synchronous generators in the system.

Definition at line 382 of file ElectricCalculation.h.

◆ m_syncMotorList

std::vector<SyncMotor*> ElectricCalculation::m_syncMotorList

protected

List of synchronous motors in the system.

Definition at line 384 of file ElectricCalculation.h.

◆ m_transformerList

std::vector<Transformer*> ElectricCalculation::m_transformerList

protected

List of transformers in the system.

Definition at line 386 of file ElectricCalculation.h.

The documentation for this class was generated from the following files:

Project/simulation/ElectricCalculation.h
Project/simulation/ElectricCalculation.cpp

Public Member Functions

Protected Member Functions

Protected Attributes

Detailed Description

Constructor & Destructor Documentation

◆ ElectricCalculation()

◆ ~ElectricCalculation()

Member Function Documentation

◆ ABCtoDQ0()

◆ CalculateEMTElementsAdmittance()

◆ CalculateEMTElementsPower()

◆ CalculateEMTPowerError()

◆ ComplexMatrixTimesVector()

◆ DistributeReactivePower()

◆ DQ0toABC()

◆ GaussianElimination() [1/2]

◆ GaussianElimination() [2/2]

◆ GetBusList()

◆ GetCapacitorList()

◆ GetElementsFromList()

◆ GetEMTElementList()

◆ GetHarmCurrentList()

◆ GetIndMotorList()

◆ GetInductorList()

◆ GetLineList()

◆ GetLoadList()

◆ GetLUDecomposition()

◆ GetMachineModel()

◆ GetNextConnection()

◆ GetParentBus() [1/2]

◆ GetParentBus() [2/2]

◆ GetPowerElementList()

◆ GetSyncGeneratorList()

◆ GetSyncMotorList()

◆ GetTransformerList()

◆ GetYBus()

◆ InvertMatrix()

◆ LUEvaluate()

◆ UpdateElementsPowerFlow()

Member Data Documentation

◆ m_busList

◆ m_capacitorList

◆ m_emtElementList

◆ m_harmCurrentList

◆ m_indMotorList

◆ m_inductorList

◆ m_lineList

◆ m_loadList

◆ m_powerElementList

◆ m_syncGeneratorList

◆ m_syncMotorList

◆ m_transformerList