process_cooling.h

#pragma once
 
#include <algorithm>
#include <cmath>
#include <fstream>
#include <stdexcept>
#include <string>
#include <utility>
#include <vector>
 
using namespace std;
 
constexpr int MONTHS             = 12;
constexpr int LOAD_NUM           = 11;
constexpr int HOURS_IN_YEAR      = 8760;
constexpr int CHILLER_COEFFS_CNT = 7; // 7 columns for each chiller, polynomial to either 3rd (for FullLoadEffKnown or custom chiller) or 5th degree
 
class ProcessCooling {
  public:
    enum RefrigerantType { R_11, R_123, R_12, R_134a, R_22, R_717 };
 
    enum ACSourceLocation { Inside, Outside };
 
    enum CoolingSystemType { Water, Air };
 
    enum CellFanType { AxialFan, CentrifugalFan };
 
    enum TowerSizedBy { Tonnage, Fan_HP, Unknown };
 
    enum ChillerCompressorType { Centrifugal, Screw, Reciprocating };
 
    enum FanMotorSpeedType { One, Two, Variable };
 
    struct ChillerOutput {
        ChillerOutput(vector<vector<double>> efficiency, vector<vector<double>> hours, vector<vector<double>> power,
                      vector<vector<double>> energy)
            : efficiency(std::move(efficiency)), hours(std::move(hours)), power(std::move(power)),
              energy(std::move(energy)) {}
 
        vector<vector<double>> efficiency;
        vector<vector<double>> hours;
        vector<vector<double>> power;
        vector<vector<double>> energy;
    };
 
    struct ChillerPumpingEnergyOutput {
        explicit ChillerPumpingEnergyOutput(vector<double> pumpingEnergy)
            : chillerPumpingEnergy(std::move(pumpingEnergy)) {}
 
        vector<double> chillerPumpingEnergy;
    };
 
    struct TowerOutput {
        TowerOutput() = default;
 
        TowerOutput(vector<double> hours, vector<double> energy) : hours(std::move(hours)), energy(std::move(energy)) {}
 
        vector<double> tempBins = {35, 45, 55, 65, 75, 75};
        vector<double> hours;
        vector<double> energy;
    };
 
    struct WaterCooledSystemInput {
        WaterCooledSystemInput() = default;
 
        WaterCooledSystemInput(double CHWT, bool useFreeCooling, double HEXApproachTemp, bool constantCWT, double CWT,
                               bool CWVariableFlow, double CWFlowRate, double CWTFollow)
            : CHWT(CHWT), useFreeCooling(useFreeCooling), HEXApproachTemp(HEXApproachTemp), constantCWT(constantCWT),
              CWT(CWT), CWVariableFlow(CWVariableFlow), CWFlowRate(CWFlowRate), CWTFollow(CWTFollow) {
            isWaterCooled = true;
        }
 
        double CHWT            = 44;
        bool   useFreeCooling  = false;
        double HEXApproachTemp = 0;
        bool   constantCWT     = true;
        double CWT             = 85;
        bool   CWVariableFlow  = true;
        double CWFlowRate      = 3;
        double CWTFollow       = 0;
        bool   isWaterCooled   = false;
    };
 
    struct AirCooledSystemInput {
        AirCooledSystemInput() = default;
 
        AirCooledSystemInput(double CHWT, double OADT, ACSourceLocation ACSource, double indoorTemp, double CWTFollow)
            : CHWT(CHWT), OADT(OADT), ACSource(ACSource), indoorTemp(indoorTemp), CWTFollow(CWTFollow) {
            isAirCooled = true;
        }
 
        double           CHWT        = 44;
        double           OADT        = 95;
        ACSourceLocation ACSource    = ACSourceLocation::Outside;
        double           indoorTemp  = 75;
        double           CWTFollow   = 0;
        bool             isAirCooled = false;
    };
 
    struct PumpInput {
        PumpInput(bool variableFlow, double flowRate, double efficiency, double motorSize, double motorEfficiency)
            : variableFlow(variableFlow), flowRate(flowRate), efficiency(efficiency * 100), motorSize(motorSize),
              motorEfficiency(motorEfficiency * 100) {}
 
        bool   variableFlow;
        double flowRate;
        double efficiency;
        double motorSize;
        double motorEfficiency;
    };
 
    struct TowerInput {
        TowerInput() = default;
 
        TowerInput(int numTower, int numFanPerTower_Cells, FanMotorSpeedType fanSpeedType, TowerSizedBy towerSizing,
                   CellFanType towerCellFanType, double cellFanHP, double tonnage)
            : numTower(numTower), numFanPerTower_Cells(numFanPerTower_Cells), fanSpeedType(fanSpeedType), towerSizing(towerSizing),
              towerCellFanType(towerCellFanType), fanHP(cellFanHP), tonnage(tonnage) {}
 
        int               numTower;
        int               numFanPerTower_Cells;
        FanMotorSpeedType fanSpeedType;
        TowerSizedBy      towerSizing;
        CellFanType       towerCellFanType;
        double            fanHP;
        double            tonnage;
    };
 
    struct ChillerInput {
        ChillerInput(ChillerCompressorType chillerType, double capacity, bool isFullLoadEffKnown, double fullLoadEff,
                     double age, bool installVSD, bool useARIMonthlyLoadSchedule, vector<vector<double>> monthlyLoads)
            : chillerType(chillerType), capacity(capacity), isFullLoadEffKnown(isFullLoadEffKnown),
              fullLoadEff(fullLoadEff), age(age), installVSD(installVSD),
              useARIMonthlyLoadSchedule(useARIMonthlyLoadSchedule), monthlyLoads(std::move(monthlyLoads)),
              isCustomChiller(false), loadAtPercent({}), kwPerTonLoads({}), changeRefrig(false),
              currentRefrig(RefrigerantType::R_11), proposedRefrig(RefrigerantType::R_11) {
            InitNonVaryingMonthlyLoad();
        }
 
        ChillerInput(ChillerCompressorType chillerType, double capacity, bool isFullLoadEffKnown, double fullLoadEff,
                     double age, bool installVSD, bool useARIMonthlyLoadSchedule, vector<vector<double>> monthlyLoads,
                     bool changeRefrig, RefrigerantType currentRefrig, RefrigerantType proposedRefrig)
            : chillerType(chillerType), capacity(capacity), isFullLoadEffKnown(isFullLoadEffKnown),
              fullLoadEff(fullLoadEff), age(age), installVSD(installVSD),
              useARIMonthlyLoadSchedule(useARIMonthlyLoadSchedule), monthlyLoads(std::move(monthlyLoads)),
              isCustomChiller(false), loadAtPercent({}), kwPerTonLoads({}), changeRefrig(changeRefrig),
              currentRefrig(currentRefrig), proposedRefrig(proposedRefrig) {
            InitNonVaryingMonthlyLoad();
        }
 
        ChillerInput(ChillerCompressorType chillerType, double capacity, bool isFullLoadEffKnown, double fullLoadEff,
                     double age, bool installVSD, bool useARIMonthlyLoadSchedule, vector<vector<double>> monthlyLoads,
                     RefrigerantType currentRefrig, RefrigerantType proposedRefrig)
            : chillerType(chillerType), capacity(capacity), isFullLoadEffKnown(isFullLoadEffKnown),
              fullLoadEff(fullLoadEff), age(age), installVSD(installVSD),
              useARIMonthlyLoadSchedule(useARIMonthlyLoadSchedule), monthlyLoads(std::move(monthlyLoads)),
              isCustomChiller(false), loadAtPercent({}), kwPerTonLoads({}), changeRefrig(true),
              currentRefrig(currentRefrig), proposedRefrig(proposedRefrig) {
            InitNonVaryingMonthlyLoad();
        }
 
        ChillerInput(ChillerCompressorType chillerType, double capacity, bool isFullLoadEffKnown, double fullLoadEff,
                     double age, bool installVSD, bool useARIMonthlyLoadSchedule, vector<vector<double>> monthlyLoads,
                     vector<double> loadAtPercent, vector<double> kwPerTonLoads)
            : chillerType(chillerType), capacity(capacity), isFullLoadEffKnown(isFullLoadEffKnown),
              fullLoadEff(fullLoadEff), age(age), installVSD(installVSD),
              useARIMonthlyLoadSchedule(useARIMonthlyLoadSchedule), monthlyLoads(std::move(monthlyLoads)),
              isCustomChiller(true), loadAtPercent(std::move(loadAtPercent)), kwPerTonLoads(std::move(kwPerTonLoads)),
              changeRefrig(false), currentRefrig(RefrigerantType::R_11), proposedRefrig(RefrigerantType::R_11) {
            InitNonVaryingMonthlyLoad();
            SetCustomCoefficient();
        }
 
        ChillerInput(ChillerCompressorType chillerType, double capacity, bool isFullLoadEffKnown, double fullLoadEff,
                     double age, bool installVSD, bool useARIMonthlyLoadSchedule, vector<vector<double>> monthlyLoads,
                     vector<double> loadAtPercent, vector<double> kwPerTonLoads, RefrigerantType currentRefrig,
                     RefrigerantType proposedRefrig)
            : chillerType(chillerType), capacity(capacity), isFullLoadEffKnown(isFullLoadEffKnown),
              fullLoadEff(fullLoadEff), age(age), installVSD(installVSD),
              useARIMonthlyLoadSchedule(useARIMonthlyLoadSchedule), monthlyLoads(std::move(monthlyLoads)),
              isCustomChiller(true), loadAtPercent(std::move(loadAtPercent)), kwPerTonLoads(std::move(kwPerTonLoads)),
              changeRefrig(true), currentRefrig(currentRefrig), proposedRefrig(proposedRefrig) {
            InitNonVaryingMonthlyLoad();
            SetCustomCoefficient();
        }
 
        ChillerCompressorType  chillerType;
        double                 capacity;
        bool                   isFullLoadEffKnown;
        double                 fullLoadEff;
        double                 age;
        bool                   installVSD;
        bool                   useARIMonthlyLoadSchedule;
        vector<vector<double>> monthlyLoads;
 
        bool           isCustomChiller;
        vector<double> loadAtPercent;
        vector<double> kwPerTonLoads;
 
        bool            changeRefrig = false;
        RefrigerantType currentRefrig;
        RefrigerantType proposedRefrig;
 
        vector<double> customCoeffs;
 
      private:
        void InitNonVaryingMonthlyLoad() {
            if (monthlyLoads.size() == 1) {
                auto monthlyLoad = monthlyLoads[0];
 
                for (int i = 0; i < 11; i++) {
                    monthlyLoads.push_back(monthlyLoad);
                }
            }
        }
 
        void SetCustomCoefficient() {
            if (loadAtPercent.empty() || kwPerTonLoads.empty() || loadAtPercent.size() != kwPerTonLoads.size()) {
                throw std::runtime_error("Invalid input provided for loadAtPercent or kwPerTonLoads, should be non empty and have same number of elements.");
            }
 
            auto size = static_cast<int>(loadAtPercent.size());
 
            // % loading in ascending order (25, 50, 75, 100)
            if (loadAtPercent[0] > loadAtPercent[size-1]) {
                std::reverse(loadAtPercent.begin(), loadAtPercent.end());
                std::reverse(kwPerTonLoads.begin(), kwPerTonLoads.end());
            }
 
            vector<double> x(size, 0);
            vector<double> y(size, 0);
            const auto kwPerTonLoadAtMaxLoad = kwPerTonLoads[size-1];
 
            if (kwPerTonLoadAtMaxLoad == 0) {
                throw std::runtime_error("% loading @ 100 % cannot be zero.");
            }
 
            for (int i = 0; i < size; i++) {
                x[i] = loadAtPercent[i]/100.0;
                y[i] = kwPerTonLoads[i] * x[i] / kwPerTonLoadAtMaxLoad;
            }
            vector<double> coeff = solveForCoefficients(x, y);
 
            size = static_cast<int>(coeff.size());
            customCoeffs.resize(size, 0);
            for (int i = 0; i < size; i++) {
                customCoeffs[i] = coeff[i];
            }
        }
 
        static vector<double> solveForCoefficients(const vector<double>& x, vector<double> y) {
            if (x.empty() || x.size() != y.size())
                return {};
 
            const int n = static_cast<int>(x.size());
 
            vector<vector<double>> a(n, vector<double>(n, 0));
            for (int i = 0; i < n; ++i) {
                for (int j = 0; j < n; ++j) {
                    a[i][j] = pow(x[i], n - 1 - j);
                }
            }
 
            // Forward Elimination
            for (int k = 0; k < n - 1; ++k) {
                for (int i = k + 1; i < n; ++i) {
                    double factor = a[i][k] / a[k][k];
                    for (int j = k + 1; j < n; ++j) {
                        a[i][j] -= factor * a[k][j];
                    }
                    y[i] -= factor * y[k];
                }
            }
 
            // Back Substitution
            vector<double> coeff(n, 0);
            coeff[n - 1] = y[n - 1] / a[n - 1][n - 1];
            for (int i = n - 2; i >= 0; --i) {
                double sum = y[i];
                for (int j = i + 1; j < n; ++j) {
                    sum -= a[i][j] * coeff[j];
                }
                coeff[i] = sum / a[i][i];
            }
 
            return coeff;
        }
    };
 
    ~ProcessCooling() = default;
 
    ProcessCooling(const vector<int>& systemOperationAnnualHours, const vector<double>& weatherDryBulbHourlyTemp,
                   const vector<double>& weatherWetBulbHourlyTemp, const vector<ChillerInput>& chillerInputList,
                   const TowerInput& towerInput, const WaterCooledSystemInput& waterCooledSystemInput)
        : ProcessCooling(systemOperationAnnualHours, weatherDryBulbHourlyTemp, weatherWetBulbHourlyTemp,
                         chillerInputList, {}, towerInput, waterCooledSystemInput) {}
 
    ProcessCooling(const vector<int>& systemOperationAnnualHours, const vector<double>& weatherDryBulbHourlyTemp,
                   const vector<double>& weatherWetBulbHourlyTemp, const vector<ChillerInput>& chillerInputList,
                   const AirCooledSystemInput& airCooledSystemInput)
        : ProcessCooling(systemOperationAnnualHours, weatherDryBulbHourlyTemp, weatherWetBulbHourlyTemp,
                         chillerInputList, airCooledSystemInput, {}, {}) {}
 
    TowerOutput calculateTowerEnergy();
 
    ChillerOutput calculateChillerEnergy();
 
    ChillerPumpingEnergyOutput calculatePumpEnergy(PumpInput pump) const;
 
    static vector<int> getSysOpAnnualHours(const vector<int>& weeklyOpStartHour, const vector<int>& weeklyOpStopHour, const vector<int>& monthlyOpMaxHour);
 
    vector<double> getChillerEfficiencyCoeffs(int chillerIndex) const;
 
    vector<double> getChillerEnergyEfficiency(int chillerIndex, const vector<double>& loadAtPercent) const;
 
  private:
    ProcessCooling(const vector<int>& systemOperationAnnualHours, const vector<double>& weatherDryBulbHourlyTemp,
                   const vector<double>& weatherWetBulbHourlyTemp, const vector<ChillerInput>& chillerInputList,
                   const AirCooledSystemInput& airCooledSystemInput, const TowerInput& towerInput,
                   const WaterCooledSystemInput& waterCooledSystemInput);
 
    vector<int>    systemOperationAnnual;
    vector<double> dryBulbHourlyTemp;
    vector<double> wetBulbHourlyTemp;
 
    TowerInput             tower {};
    WaterCooledSystemInput waterCooledSystem;
    AirCooledSystemInput   airCooledSystem;
    CoolingSystemType      coolingType;
 
    double         FCTemp = 0; // Free Cooling Temperature
    vector<double> CWTHourly;
 
    int                    numChillers;
    vector<ChillerInput>   chillers;
    vector<vector<double>> chillerHourlyLoad;
    vector<vector<double>> chillerHourlyLoadOperational;
    vector<vector<double>> chillerEfficiencyCoeffs;
    vector<vector<double>> chillerHourlyEfficiencyARI;
    vector<vector<double>> chillerHourlyEfficiency;
    vector<vector<double>> chillerHourlyPower;
 
    int getChillerCapacityIndex(ChillerCompressorType chillerType, double capacity) const;
 
    double getChillerEffAtLoad(int c, double load, bool isFullLoadEffKnown) const;
 
    void annualChillerLoadProfile();
 
    void annualChillerEfficiencyProfileARI();
 
    void annualChillerEfficiencyProfile();
 
    void annualChillerPowerProfile();
 
    void setTowerFanHPTonnage();
 
    double getPercentFanPower(double wetBulbTemp, double percentWaterFlow, double range, double desiredApproach,
                              int yearHourIndex);
 
    double getPercentWaterFlow(int yearHourIndex) const;
 
    double getRange(int yearHourIndex) const;
 
    double getApproach(double wetBulbTemp) const;
 
    double modifyPercentFanPower(double percentFanPower) const;
 
    double getWeightedAverageChillerLoad(int yearHourIndex) const;
 
    double getChillerTonnageTotal() const;
 
    static double getCubeRoot(double number);
 
    static double getPumpHP(double power);
};