PowerFlowApplication.C

Go to the documentation of this file.

/*
 * Copyright (c) 2024 Saud Zahir
 *
 * This file is part of deltaFlow.
 *
 * deltaFlow is free software; you can redistribute it and/or
 * modify it under the terms of the GNU General Public
 * License as published by the Free Software Foundation; either
 * version 3 of the License, or (at your option) any later version.
 *
 * deltaFlow is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * General Public License for more details.
 *
 * You should have received a copy of the GNU General Public
 * License along with deltaFlow.  If not, see
 * <https://www.gnu.org/licenses/>.
 */
 
#include "PowerFlowApplication.H"
 
#include <chrono>
#include <cmath>
#include <complex>
#include <cstdlib>
#include <utility>
 
#include "Admittance.H"
#include "Argparse.H"
#include "DatFileWriter.H"
#include "Display.H"
#include "Logger.H"
#include "LoggingSolverObserver.H"
#include "MsgFileWriter.H"
#include "OutFileWriter.H"
#include "Qlim.H"
#include "ReaderFactory.H"
#include "SolverFactory.H"
#include "StaFileWriter.H"
#include "Writer.H"
 
PowerFlowApplication::PowerFlowApplication(PowerFlowOptions options)
    : options_(std::move(options)) {}
 
void PowerFlowApplication::loadInput() {
    std::string inputFile = options_.inputFile;
    InputFormat format = options_.inputFormat;
 
    reader_ = ReaderFactory::create(format);
    formatName_ = reader_->name();
 
    LOG_INFO("Reading {} file: {}", formatName_, inputFile);
 
    reader_->read(inputFile);
 
    busData_ = reader_->getBusData();
    branchData_ = reader_->getBranchData();
 
    if (busData_.ID.size() == 0 || branchData_.From.size() == 0) {
        LOG_ERROR("No bus or branch data found in '{}'. Check the file exists and is valid.", inputFile);
        std::exit(1);
    }
 
    N_ = busData_.ID.size();
    nBranch_ = branchData_.From.size();
 
    LOG_INFO("Model: {} buses, {} branches", N_, nBranch_);
 
    int nSlack = 0, nPV = 0, nPQ = 0;
    for (int i = 0; i < N_; ++i) {
        if (busData_.Type(i) == 1) nSlack++;
        else if (busData_.Type(i) == 2) nPV++;
        else nPQ++;
    }
    LOG_DEBUG("Bus types: {} Slack, {} PV, {} PQ", nSlack, nPV, nPQ);
}
 
void PowerFlowApplication::prepareModel() {
    LoggingSolverObserver observer;
    Y_ = Admittance::build(busData_, branchData_, &observer);
    LOG_DEBUG("Admittance matrix computed ({}x{})", N_, N_);
 
    V_.resize(N_);
    delta_ = Eigen::VectorXd::Zero(N_);
 
    // Flat start: PQ buses -> V=1.0, delta=0 for all; PV/Slack keep file voltage
    for (int i = 0; i < N_; ++i) {
        if (busData_.Type(i) == 3) {  // PQ bus
            V_(i) = 1.0;
        } else {
            V_(i) = busData_.V(i);  // PV/Slack keep specified voltage
        }
    }
 
    typeBus_ = busData_.Type;
}
 
void PowerFlowApplication::solve() {
    SolverType solverType = options_.solverType;
 
    auto solverImpl = SolverFactory::create(solverType);
    solverName_ = solverImpl->name();
 
    Eigen::MatrixXd G = Y_.array().real().matrix();
    Eigen::MatrixXd B = Y_.array().imag().matrix();
 
    int maxIter = options_.maxIterations;
    double tolerance = options_.tolerance;
 
    LOG_INFO("Starting {} solver ...", solverName_);
 
    LoggingSolverObserver observer;
    bool Q_lim_status = true;
 
    while (Q_lim_status) {
        Eigen::VectorXd Ps = busData_.Pg - busData_.Pl;
        Eigen::VectorXd Qs = busData_.Qg - busData_.Ql;
 
        std::vector<int> pq_indices;
        std::vector<int> pv_indices;
 
        for (int i = 0; i < N_; ++i) {
            if (typeBus_(i) == 3) pq_indices.push_back(i);
            else if (typeBus_(i) == 2) pv_indices.push_back(i);
        }
 
        SolverContext ctx;
        ctx.Y = &Y_;
        ctx.G = &G;
        ctx.B = &B;
        ctx.V = &V_;
        ctx.delta = &delta_;
        ctx.type_bus = &typeBus_;
        ctx.pq_bus_id = &pq_indices;
        ctx.P = &Ps;
        ctx.Q = &Qs;
        ctx.N = N_;
        ctx.n_pq = static_cast<int>(pq_indices.size());
        ctx.maxIter = maxIter;
        ctx.tolerance = tolerance;
        ctx.omega = options_.relaxation;
        ctx.iterHistory = &iterationHistory_;
        ctx.observer = &observer;
 
        bool converged = solverImpl->solve(ctx);
        finalConverged_ = converged;
 
        if (!converged) {
            LOG_ERROR("{} solver failed to converge.", solverName_);
            break;
        }
 
        Q_lim_status = Qlim::check(V_, delta_, typeBus_, G, B, busData_, pv_indices, N_, &observer);
 
        if (Q_lim_status) {
            LOG_DEBUG("Re-running {} with updated bus types ...", solverName_);
        }
    }
 
    if (!iterationHistory_.empty()) {
        totalIterations_ = iterationHistory_.back().first;
        finalError_ = iterationHistory_.back().second;
    }
}
 
void PowerFlowApplication::finalizeResults() {
    Eigen::VectorXcd Vc(N_);
    for (int i = 0; i < N_; ++i)
        Vc(i) = std::polar(V_(i), delta_(i));
 
    Eigen::VectorXd P_net = busData_.Pg - busData_.Pl;
    Eigen::VectorXd Q_net = busData_.Qg - busData_.Ql;
 
    // Recalculate slack bus power injection
    for (int i = 0; i < N_; ++i) {
        if (busData_.Type(i) == 1) {  // Slack
            std::complex<double> Ii = Y_.row(i) * Vc;
            std::complex<double> Si = Vc(i) * std::conj(Ii);
            P_net(i) = Si.real();
            Q_net(i) = Si.imag();
        }
    }
 
    // Recalculate reactive power for PV buses
    for (int i = 0; i < N_; ++i) {
        if (busData_.Type(i) == 2) {  // PV
            std::complex<double> Ii = Y_.row(i) * Vc;
            Q_net(i) = -std::imag(std::conj(Vc(i)) * Ii);
        }
    }
 
    for (int i = 0; i < N_; ++i) {
        busData_.V(i) = std::abs(Vc(i));
        busData_.delta(i) = std::arg(Vc(i)) * 180.0 / M_PI;
        busData_.Pg(i) = P_net(i) + busData_.Pl(i);
        busData_.Qg(i) = Q_net(i) + busData_.Ql(i);
    }
 
    double PLoss = busData_.Pg.sum() - busData_.Pl.sum();
    double QLoss = busData_.Qg.sum() - busData_.Ql.sum();
    LOG_DEBUG("Total real power loss: {:.6f} p.u.", PLoss);
    LOG_DEBUG("Total reactive power loss: {:.6f} p.u.", QLoss);
}
 
void PowerFlowApplication::writeFailureStatus(double elapsedSec) const {
    StaFileWriter staWriter;
    staWriter.write(options_.jobName, options_.inputFile, solverName_, formatName_,
        N_, nBranch_, options_.maxIterations, 0.0, options_.tolerance, false, elapsedSec);
}
 
void PowerFlowApplication::writeOutputs(double elapsedSec) const {
    dispBusData(busData_);
    dispLineFlow(busData_, branchData_, Y_);
 
    writeOutputCSV(busData_);
 
    std::string jobName = options_.jobName;
    std::string inputFile = options_.inputFile;
    double tolerance = options_.tolerance;
 
    OutFileWriter outWriter;
    outWriter.write(jobName, inputFile, solverName_, formatName_,
        busData_, branchData_, Y_, totalIterations_, finalError_, tolerance, elapsedSec);
 
    StaFileWriter staWriter;
    staWriter.write(jobName, inputFile, solverName_, formatName_,
        N_, nBranch_, totalIterations_, finalError_, tolerance, finalConverged_, elapsedSec);
 
    DatFileWriter datWriter;
    datWriter.write(jobName, inputFile, solverName_, formatName_,
        busData_, branchData_, iterationHistory_, totalIterations_, finalError_,
        tolerance, finalConverged_, elapsedSec);
 
    MsgFileWriter msgWriter;
    msgWriter.write(jobName, solverName_, iterationHistory_, tolerance, finalConverged_);
}
 
int PowerFlowApplication::run() {
    auto startTime = std::chrono::high_resolution_clock::now();
 
    LOG_DEBUG("Job name     :: {}", options_.jobName);
    LOG_DEBUG("Input file   :: {}", options_.inputFile);
    LOG_DEBUG("Tolerance    :: {:.6e}", options_.tolerance);
    LOG_DEBUG("Max iter     :: {}", options_.maxIterations);
 
    loadInput();
    prepareModel();
    solve();
 
    if (!finalConverged_) {
        auto endTime = std::chrono::high_resolution_clock::now();
        double elapsedSec = std::chrono::duration<double>(endTime - startTime).count();
        writeFailureStatus(elapsedSec);
        return 1;
    }
 
    finalizeResults();
 
    auto endTime = std::chrono::high_resolution_clock::now();
    double elapsedSec = std::chrono::duration<double>(endTime - startTime).count();
 
    writeOutputs(elapsedSec);
 
    fmt::print("\n");
    fmt::print(fg(Display::LOGO_COLOR) | fmt::emphasis::bold,
        "   THE ANALYSIS HAS BEEN COMPLETED SUCCESSFULLY\n");
    fmt::print("\n");
    fmt::print("   Elapsed time : {:.3f} sec\n", elapsedSec);
    fmt::print("\n");
 
    return 0;
}