EMT_Ph3_SSN_Capacitor.cpp

/* Copyright 2017-2021 Institute for Automation of Complex Power Systems,
 *                     EONERC, RWTH Aachen University
 *
 * This Source Code Form is subject to the terms of the Mozilla Public
 * License, v. 2.0. If a copy of the MPL was not distributed with this
 * file, You can obtain one at https://mozilla.org/MPL/2.0/.
 *********************************************************************************/
 
#include <dpsim-models/EMT/EMT_Ph3_SSN_Capacitor.h>
 
using namespace CPS;
 
EMT::Ph3::SSN::Capacitor::Capacitor(String uid, String name,
                                    Logger::Level logLevel)
    : MNASimPowerComp<Real>(uid, name, true, true, logLevel),
      Base::Ph3::Capacitor(mAttributes) {
  mPhaseType = PhaseType::ABC;
  setVirtualNodeNumber(1);
  setTerminalNumber(2);
  mHistoricVoltage = Matrix::Zero(3, 1);
  **mIntfVoltage = Matrix::Zero(3, 1);
  **mIntfCurrent = Matrix::Zero(3, 1);
}
 
SimPowerComp<Real>::Ptr EMT::Ph3::SSN::Capacitor::clone(String name) {
  auto copy = Capacitor::make(name, mLogLevel);
  copy->setParameters(**mCapacitance);
  return copy;
}
void EMT::Ph3::SSN::Capacitor::initializeFromNodesAndTerminals(Real frequency) {
 
  Real omega = 2 * PI * frequency;
  MatrixComp admittance = MatrixComp::Zero(3, 3);
  admittance << Complex(0, omega * (**mCapacitance)(0, 0)),
      Complex(0, omega * (**mCapacitance)(0, 1)),
      Complex(0, omega * (**mCapacitance)(0, 2)),
      Complex(0, omega * (**mCapacitance)(1, 0)),
      Complex(0, omega * (**mCapacitance)(1, 1)),
      Complex(0, omega * (**mCapacitance)(1, 2)),
      Complex(0, omega * (**mCapacitance)(2, 0)),
      Complex(0, omega * (**mCapacitance)(2, 1)),
      Complex(0, omega * (**mCapacitance)(2, 2));
 
  MatrixComp vInitABC = Matrix::Zero(3, 1);
  vInitABC(0, 0) = RMS3PH_TO_PEAK1PH * initialSingleVoltage(1) -
                   RMS3PH_TO_PEAK1PH * initialSingleVoltage(0);
  vInitABC(1, 0) = vInitABC(0, 0) * SHIFT_TO_PHASE_B;
  vInitABC(2, 0) = vInitABC(0, 0) * SHIFT_TO_PHASE_C;
  **mIntfVoltage = vInitABC.real();
  **mIntfCurrent = (admittance * vInitABC).real();
 
  SPDLOG_LOGGER_INFO(mSLog,
                     "\nCapacitance [F]: {:s}"
                     "\nAdmittance [S]: {:s}",
                     Logger::matrixToString(**mCapacitance),
                     Logger::matrixCompToString(admittance));
  SPDLOG_LOGGER_INFO(
      mSLog,
      "\n--- Initialization from powerflow ---"
      "\nVoltage across: {:s}"
      "\nCurrent: {:s}"
      "\nTerminal 0 voltage: {:s}"
      "\nTerminal 1 voltage: {:s}"
      "\n--- Initialization from powerflow finished ---",
      Logger::matrixToString(**mIntfVoltage),
      Logger::matrixToString(**mIntfCurrent),
      Logger::phasorToString(RMS3PH_TO_PEAK1PH * initialSingleVoltage(0)),
      Logger::phasorToString(RMS3PH_TO_PEAK1PH * initialSingleVoltage(1)));
}
 
void EMT::Ph3::SSN::Capacitor::mnaCompInitialize(
    Real omega, Real timeStep, Attribute<Matrix>::Ptr leftVector) {
  updateMatrixNodeIndices();
 
  mHistoricVoltage = mDufourBKHat * **mIntfCurrent + **mIntfVoltage;
  mDufourBKHat = timeStep * (2.0 * **mCapacitance).inverse();
  mDufourWKN = mDufourBKHat;
 
  **mRightVector = Matrix::Zero(leftVector->get().rows(), 1);
}
 
void EMT::Ph3::SSN::Capacitor::mnaCompApplySystemMatrixStamp(
    SparseMatrixRow &systemMatrix) {
  if (terminalNotGrounded(0)) {
    Math::addToMatrixElement(systemMatrix, matrixNodeIndex(0, 0),
                             mVirtualNodes[0]->matrixNodeIndex(PhaseType::A),
                             -1);
    Math::addToMatrixElement(systemMatrix,
                             mVirtualNodes[0]->matrixNodeIndex(PhaseType::A),
                             matrixNodeIndex(0, 0), -1);
 
    Math::addToMatrixElement(systemMatrix, matrixNodeIndex(0, 1),
                             mVirtualNodes[0]->matrixNodeIndex(PhaseType::B),
                             -1);
    Math::addToMatrixElement(systemMatrix,
                             mVirtualNodes[0]->matrixNodeIndex(PhaseType::B),
                             matrixNodeIndex(0, 1), -1);
 
    Math::addToMatrixElement(systemMatrix, matrixNodeIndex(0, 2),
                             mVirtualNodes[0]->matrixNodeIndex(PhaseType::C),
                             -1);
    Math::addToMatrixElement(systemMatrix,
                             mVirtualNodes[0]->matrixNodeIndex(PhaseType::C),
                             matrixNodeIndex(0, 2), -1);
  }
  if (terminalNotGrounded(1)) {
    Math::addToMatrixElement(systemMatrix, matrixNodeIndex(1, 0),
                             mVirtualNodes[0]->matrixNodeIndex(PhaseType::A),
                             1);
    Math::addToMatrixElement(systemMatrix,
                             mVirtualNodes[0]->matrixNodeIndex(PhaseType::A),
                             matrixNodeIndex(1, 0), 1);
 
    Math::addToMatrixElement(systemMatrix, matrixNodeIndex(1, 1),
                             mVirtualNodes[0]->matrixNodeIndex(PhaseType::B),
                             1);
    Math::addToMatrixElement(systemMatrix,
                             mVirtualNodes[0]->matrixNodeIndex(PhaseType::B),
                             matrixNodeIndex(1, 1), 1);
 
    Math::addToMatrixElement(systemMatrix, matrixNodeIndex(1, 2),
                             mVirtualNodes[0]->matrixNodeIndex(PhaseType::C),
                             1);
    Math::addToMatrixElement(systemMatrix,
                             mVirtualNodes[0]->matrixNodeIndex(PhaseType::C),
                             matrixNodeIndex(1, 2), 1);
  }
  //mesh equations are independent from grounded terminals
  Math::addToMatrixElement(
      systemMatrix, mVirtualNodes[0]->matrixNodeIndex(PhaseType::A),
      mVirtualNodes[0]->matrixNodeIndex(PhaseType::A), -mDufourWKN(0, 0));
  Math::addToMatrixElement(
      systemMatrix, mVirtualNodes[0]->matrixNodeIndex(PhaseType::A),
      mVirtualNodes[0]->matrixNodeIndex(PhaseType::B), -mDufourWKN(0, 1));
  Math::addToMatrixElement(
      systemMatrix, mVirtualNodes[0]->matrixNodeIndex(PhaseType::A),
      mVirtualNodes[0]->matrixNodeIndex(PhaseType::C), -mDufourWKN(0, 2));
 
  Math::addToMatrixElement(
      systemMatrix, mVirtualNodes[0]->matrixNodeIndex(PhaseType::B),
      mVirtualNodes[0]->matrixNodeIndex(PhaseType::A), -mDufourWKN(1, 0));
  Math::addToMatrixElement(
      systemMatrix, mVirtualNodes[0]->matrixNodeIndex(PhaseType::B),
      mVirtualNodes[0]->matrixNodeIndex(PhaseType::B), -mDufourWKN(1, 1));
  Math::addToMatrixElement(
      systemMatrix, mVirtualNodes[0]->matrixNodeIndex(PhaseType::B),
      mVirtualNodes[0]->matrixNodeIndex(PhaseType::C), -mDufourWKN(1, 2));
 
  Math::addToMatrixElement(
      systemMatrix, mVirtualNodes[0]->matrixNodeIndex(PhaseType::C),
      mVirtualNodes[0]->matrixNodeIndex(PhaseType::A), -mDufourWKN(2, 0));
  Math::addToMatrixElement(
      systemMatrix, mVirtualNodes[0]->matrixNodeIndex(PhaseType::C),
      mVirtualNodes[0]->matrixNodeIndex(PhaseType::B), -mDufourWKN(2, 1));
  Math::addToMatrixElement(
      systemMatrix, mVirtualNodes[0]->matrixNodeIndex(PhaseType::C),
      mVirtualNodes[0]->matrixNodeIndex(PhaseType::C), -mDufourWKN(2, 2));
}
 
void EMT::Ph3::SSN::Capacitor::mnaCompApplyRightSideVectorStamp(
    Matrix &rightVector) {
  mHistoricVoltage = mDufourBKHat * **mIntfCurrent + **mIntfVoltage;
  Math::setVectorElement(rightVector,
                         mVirtualNodes[0]->matrixNodeIndex(PhaseType::A),
                         mHistoricVoltage(0, 0));
  Math::setVectorElement(rightVector,
                         mVirtualNodes[0]->matrixNodeIndex(PhaseType::B),
                         mHistoricVoltage(1, 0));
  Math::setVectorElement(rightVector,
                         mVirtualNodes[0]->matrixNodeIndex(PhaseType::C),
                         mHistoricVoltage(2, 0));
}
 
void EMT::Ph3::SSN::Capacitor::mnaCompAddPreStepDependencies(
    AttributeBase::List &prevStepDependencies,
    AttributeBase::List &attributeDependencies,
    AttributeBase::List &modifiedAttributes) {
  // actually depends on C, but then we'd have to modify the system matrix anyway
  prevStepDependencies.push_back(mIntfCurrent);
  prevStepDependencies.push_back(mIntfVoltage);
  modifiedAttributes.push_back(mRightVector);
}
 
void EMT::Ph3::SSN::Capacitor::mnaCompPreStep(Real time, Int timeStepCount) {
  mnaCompApplyRightSideVectorStamp(**mRightVector);
}
 
void EMT::Ph3::SSN::Capacitor::mnaCompAddPostStepDependencies(
    AttributeBase::List &prevStepDependencies,
    AttributeBase::List &attributeDependencies,
    AttributeBase::List &modifiedAttributes,
    Attribute<Matrix>::Ptr &leftVector) {
  attributeDependencies.push_back(leftVector);
  modifiedAttributes.push_back(mIntfVoltage);
  modifiedAttributes.push_back(mIntfCurrent);
}
 
void EMT::Ph3::SSN::Capacitor::mnaCompPostStep(
    Real time, Int timeStepCount, Attribute<Matrix>::Ptr &leftVector) {
  mnaCompUpdateVoltage(**leftVector);
  mnaCompUpdateCurrent(**leftVector);
}
 
void EMT::Ph3::SSN::Capacitor::mnaCompUpdateVoltage(const Matrix &leftVector) {
  // v1 - v0
  **mIntfVoltage = Matrix::Zero(3, 1);
  if (terminalNotGrounded(1)) {
    (**mIntfVoltage)(0, 0) =
        Math::realFromVectorElement(leftVector, matrixNodeIndex(1, 0));
    (**mIntfVoltage)(1, 0) =
        Math::realFromVectorElement(leftVector, matrixNodeIndex(1, 1));
    (**mIntfVoltage)(2, 0) =
        Math::realFromVectorElement(leftVector, matrixNodeIndex(1, 2));
  }
  if (terminalNotGrounded(0)) {
    (**mIntfVoltage)(0, 0) =
        (**mIntfVoltage)(0, 0) -
        Math::realFromVectorElement(leftVector, matrixNodeIndex(0, 0));
    (**mIntfVoltage)(1, 0) =
        (**mIntfVoltage)(1, 0) -
        Math::realFromVectorElement(leftVector, matrixNodeIndex(0, 1));
    (**mIntfVoltage)(2, 0) =
        (**mIntfVoltage)(2, 0) -
        Math::realFromVectorElement(leftVector, matrixNodeIndex(0, 2));
  }
}
 
void EMT::Ph3::SSN::Capacitor::mnaCompUpdateCurrent(const Matrix &leftVector) {
  (**mIntfCurrent)(0, 0) = Math::realFromVectorElement(
      leftVector, mVirtualNodes[0]->matrixNodeIndex(PhaseType::A));
  (**mIntfCurrent)(1, 0) = Math::realFromVectorElement(
      leftVector, mVirtualNodes[0]->matrixNodeIndex(PhaseType::B));
  (**mIntfCurrent)(2, 0) = Math::realFromVectorElement(
      leftVector, mVirtualNodes[0]->matrixNodeIndex(PhaseType::C));
 
  SPDLOG_LOGGER_DEBUG(mSLog, "\nCurrent: {:s}",
                      Logger::matrixToString(**mIntfCurrent));
}