process_cooling.h

#pragma once
 
#include <cmath>
#include <fstream>
#include <iostream>
#include <sstream>
#include <string>
#include <utility>
#include <vector>
 
using namespace std;
 
const int MONTHS        = 12;
const int LOAD_NUM      = 11;
const int HOURS_IN_YEAR = 8760;
 
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 };
 
    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),
              fanHP(cellFanHP), tonnage(tonnage) {
            fanHP = getFanHP(tonnage, towerSizing, numFanPerTower_Cells, towerCellFanType, fanHP);
        }
 
        int               numTower;
        int               numFanPerTower_Cells;
        FanMotorSpeedType fanSpeedType;
        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() {
            auto           size = static_cast<int>(loadAtPercent.size());
            vector<double> x(size, 0);
            vector<double> y(size, 0);
            for (int i = 0; i < size; i++) {
                x[i] = loadAtPercent[i];
                y[i] = kwPerTonLoads[i] * loadAtPercent[i] / kwPerTonLoads[0];
            }
            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(vector<double> x, vector<double> y) {
            if (x.empty() || x.size() != y.size())
                return {};
 
            const int n = (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,
                   TowerInput towerInput, 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,
                   AirCooledSystemInput airCooledSystemInput)
        : ProcessCooling(systemOperationAnnualHours, weatherDryBulbHourlyTemp, weatherWetBulbHourlyTemp,
                         chillerInputList, airCooledSystemInput, {}, {}) {}
 
    TowerOutput calculateTowerEnergy();
 
    ChillerOutput calculateChillerEnergy();
 
    ChillerPumpingEnergyOutput calculatePumpEnergy(PumpInput pump);
 
    static vector<int> getSysOpAnnualHours(const vector<int>& weeklyOpStartHour, const vector<int>& weeklyOpStopHour, const vector<int>& monthlyOpMaxHour);
 
  private:
    ProcessCooling(const vector<int>& systemOperationAnnualHours, const vector<double>& weatherDryBulbHourlyTemp,
                   const vector<double>& weatherWetBulbHourlyTemp, const vector<ChillerInput>& chillerInputList,
                   AirCooledSystemInput airCooledSystemInput, TowerInput towerInput,
                   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>> chillerHourlyEfficiencyARI;
    vector<vector<double>> chillerHourlyEfficiency;
    vector<vector<double>> chillerHourlyPower;
 
    void annualChillerLoadProfile();
 
    void annualChillerEfficiencyProfileARI();
 
    void annualChillerEfficiencyProfile();
 
    void annualChillerPowerProfile();
 
    static double getFanHP(double tonnage, TowerSizedBy towerSizing, int fanNum, CellFanType fanType, double fanHP);
 
    double getPercentFanPower(double wetBulbTemp, double percentWaterFlow, double range, double desiredApproach,
                              int yearHourIndex);
 
    double getPercentWaterFlow(int yearHourIndex);
 
    double getRange(int yearHourIndex);
 
    double getApproach(double wetBulbTemp, double minToChillersTemp) const;
 
    double modifyPercentFanPower(double percentFanPower) const;
 
    double getWeightedAverageChillerLoad(int yearHourIndex);
 
    double getChillerTonnageTotal();
 
    static double getCubeRoot(double number);
 
    static double getPumpHP(double power);
};