SolveRoot with local Optimization

src/EnergyPlus/Autosizing/WaterHeatingCoilUASizing.cc

@@ -59,9 +59,6 @@ Real64 WaterHeatingCoilUASizer::size(EnergyPlusData &state, Real64 _originalValu
     if (!this->checkInitialized(state, errorsFound)) {
         return 0.0;
     }
-    Real64 constexpr Acc(0.0001); // Accuracy of result
-    int constexpr MaxIte(500);    // Maximum number of iterations
-    int SolFla = 0;               // Flag of solver
 
     this->preSize(state, _originalValue);
     if (this->curZoneEqNum > 0) {
@@ -73,15 +70,21 @@ Real64 WaterHeatingCoilUASizer::size(EnergyPlusData &state, Real64 _originalValu
                 Real64 UA1 = this->dataCapacityUsedForSizing;
                 // Invert the simple heating coil model: given the design inlet conditions and the design load,
                 // find the design UA.
+
                 auto f = [&state, this](Real64 const UA) {
                     state.dataWaterCoils->WaterCoil(this->dataCoilNum).UACoilVariable = UA;
                     WaterCoils::CalcSimpleHeatingCoil(state, this->dataCoilNum, this->dataFanOp, 1.0, state.dataWaterCoils->SimCalc);
                     state.dataSize->DataDesignCoilCapacity = state.dataWaterCoils->WaterCoil(this->dataCoilNum).TotWaterHeatingCoilRate;
                     return (dataCapacityUsedForSizing - state.dataWaterCoils->WaterCoil(this->dataCoilNum).TotWaterHeatingCoilRate) /
                            dataCapacityUsedForSizing;
                 };
-                General::SolveRoot(state, Acc, MaxIte, SolFla, this->autoSizedValue, f, UA0, UA1);
-                if (SolFla == -1) {
+
+                constexpr Real64 Acc = 0.0001; // Accuracy of result
+                int SolFla;
+                // Don't use SolveRoot2 (optimizer) for sizing
+                General::SolveRoot(state, Acc, 500, SolFla, this->autoSizedValue, f, UA0, UA1);
+
+                if (SolFla == General::SOLVEROOT_ERROR_ITER) {
                     errorsFound = true;
                     std::string msg = "Autosizing of heating coil UA failed for Coil:Heating:Water \"" + this->compName + "\"";
                     this->addErrorMessage(msg);
@@ -142,7 +145,8 @@ Real64 WaterHeatingCoilUASizer::size(EnergyPlusData &state, Real64 _originalValu
                         ShowContinueError(state, msg);
                     }
                     this->dataErrorsFound = true;
-                } else if (SolFla == -2) {
+
+                } else if (SolFla == General::SOLVEROOT_ERROR_INIT) {
                     this->errorType = AutoSizingResultType::ErrorType1;
                     errorsFound = true;
                     std::string msg = "Autosizing of heating coil UA failed for Coil:Heating:Water \"" + this->compName + "\"";
@@ -256,8 +260,13 @@ Real64 WaterHeatingCoilUASizer::size(EnergyPlusData &state, Real64 _originalValu
                     return (dataCapacityUsedForSizing - state.dataWaterCoils->WaterCoil(this->dataCoilNum).TotWaterHeatingCoilRate) /
                            dataCapacityUsedForSizing;
                 };
-                General::SolveRoot(state, Acc, MaxIte, SolFla, this->autoSizedValue, f, UA0, UA1);
-                if (SolFla == -1) {
+
+                constexpr Real64 Acc = 0.0001; // Necessary?
+                int SolFla;
+                // Don't use SolveRoot2 (optimizer) for sizing
+                General::SolveRoot(state, Acc, 500, SolFla, this->autoSizedValue, f, UA0, UA1);
+
+                if (SolFla == General::SOLVEROOT_ERROR_ITER) {
                     errorsFound = true;
                     std::string msg = "Autosizing of heating coil UA failed for Coil:Heating:Water \"" + this->compName + "\"";
                     this->addErrorMessage(msg);
@@ -302,7 +311,7 @@ Real64 WaterHeatingCoilUASizer::size(EnergyPlusData &state, Real64 _originalValu
                         ShowContinueError(state, msg);
                     }
                     this->dataErrorsFound = true;
-                } else if (SolFla == -2) {
+                } else if (SolFla == General::SOLVEROOT_ERROR_INIT) {
                     this->errorType = AutoSizingResultType::ErrorType1;
                     errorsFound = true;
                     std::string msg = "Autosizing of heating coil UA failed for Coil:Heating:Water \"" + this->compName + "\"";

src/EnergyPlus/DXCoils.cc

@@ -14839,6 +14839,7 @@ void CalcTwoSpeedDXCoilStandardRating(EnergyPlusData &state, int const DXCoilNum
                 Real64 OutletAirTemp = state.dataDXCoils->DXCoilOutletTemp(DXCoilNum);
                 return TempDryBulb_Leaving_Apoint - OutletAirTemp;
             };
+
         General::SolveRoot(state,
                            AccuracyTolerance,
                            MaximumIterations,

src/EnergyPlus/ExtendedHeatIndex.cc

@@ -542,8 +542,8 @@ namespace ExtendedHI {
         // RH: relative humidity in range of 0.0 to 1.0
         // The function computes the extended heat index, in Kelvinn
 
-        auto const HVACSystemRootSolverMethodBackup = state.dataRootFinder->HVACSystemRootFinding.HVACSystemRootSolverMethod;
-        state.dataRootFinder->HVACSystemRootFinding.HVACSystemRootSolverMethod = HVACSystemRootSolverAlgorithm::ShortBisectionThenRegulaFalsi;
+        RootAlgo rootAlgoBackup = state.dataRootFinder->rootAlgo;
+        state.dataRootFinder->rootAlgo = RootAlgo::ShortBisectionThenRegulaFalsi;
         EqvarName eqvar_name = EqvarName::Invalid;
         Real64 const eqvar_value = find_eqvar_name_and_value(state, Ta, RH, eqvar_name);
 
@@ -553,7 +553,7 @@ namespace ExtendedHI {
             T = 0.0;
         }
 
-        state.dataRootFinder->HVACSystemRootFinding.HVACSystemRootSolverMethod = HVACSystemRootSolverMethodBackup;
+        state.dataRootFinder->rootAlgo = rootAlgoBackup;
         return T;
     }
 

src/EnergyPlus/General.cc

@@ -193,7 +193,7 @@ void SolveRoot(const EnergyPlusData &state,
     Real64 Y1 = f(X1); // f at X1
     // check initial values
     if (Y0 * Y1 > 0) {
-        Flag = -2;
+        Flag = SOLVEROOT_ERROR_INIT;
         XRes = X0;
         return;
     }
@@ -209,43 +209,44 @@ void SolveRoot(const EnergyPlusData &state,
             break;
         }
         // new estimation
-        switch (state.dataRootFinder->HVACSystemRootFinding.HVACSystemRootSolverMethod) {
-        case HVACSystemRootSolverAlgorithm::RegulaFalsi: {
+        switch (state.dataRootFinder->rootAlgo) {
+        case RootAlgo::RegulaFalsi: {
             XTemp = (Y0 * X1 - Y1 * X0) / DY;
             break;
         }
-        case HVACSystemRootSolverAlgorithm::Bisection: {
+        case RootAlgo::Bisection: {
             XTemp = (X1 + X0) / 2.0;
             break;
         }
-        case HVACSystemRootSolverAlgorithm::RegulaFalsiThenBisection: {
-            if (NIte > state.dataRootFinder->HVACSystemRootFinding.NumOfIter) {
+        case RootAlgo::RegulaFalsiThenBisection: {
+            if (NIte > state.dataRootFinder->NumOfIter) {
                 XTemp = (X1 + X0) / 2.0;
             } else {
                 XTemp = (Y0 * X1 - Y1 * X0) / DY;
             }
             break;
         }
-        case HVACSystemRootSolverAlgorithm::BisectionThenRegulaFalsi: {
-            if (NIte <= state.dataRootFinder->HVACSystemRootFinding.NumOfIter) {
+        case RootAlgo::BisectionThenRegulaFalsi: {
+            if (NIte <= state.dataRootFinder->NumOfIter) {
                 XTemp = (X1 + X0) / 2.0;
             } else {
                 XTemp = (Y0 * X1 - Y1 * X0) / DY;
             }
             break;
         }
-        case HVACSystemRootSolverAlgorithm::Alternation: {
-            if (AltIte > state.dataRootFinder->HVACSystemRootFinding.NumOfIter) {
+        case RootAlgo::Alternation: {
+            if (AltIte > state.dataRootFinder->NumOfIter) {
                 XTemp = (X1 + X0) / 2.0;
-                if (AltIte >= 2 * state.dataRootFinder->HVACSystemRootFinding.NumOfIter) {
+
+                if (AltIte >= 2 * state.dataRootFinder->NumOfIter) {
                     AltIte = 0;
                 }
             } else {
                 XTemp = (Y0 * X1 - Y1 * X0) / DY;
             }
             break;
         }
-        case HVACSystemRootSolverAlgorithm::ShortBisectionThenRegulaFalsi: {
+        case RootAlgo::ShortBisectionThenRegulaFalsi: {
             if (NIte < 3) {
                 XTemp = (X1 + X0) / 2.0;
             } else {
@@ -270,6 +271,15 @@ void SolveRoot(const EnergyPlusData &state,
             return;
         };
 
+#ifdef GET_OUT
+        if (NIte > 20) {
+            assert(false);
+            Flag = NIte;
+            XRes = XTemp;
+            return;
+        }
+#endif // GET_OUT
+
         // OK, so we didn't converge, lets check max iterations to see if we should break early
         if (NIte > MaxIte) {
             break;
@@ -297,10 +307,70 @@ void SolveRoot(const EnergyPlusData &state,
     } // Cont
 
     // if we make it here we haven't converged, so just set the flag and leave
-    Flag = -1;
+    Flag = SOLVEROOT_ERROR_ITER;
     XRes = XTemp;
 }
 
+// A second version that does not require a payload -- use lambdas
+Real64 SolveRoot2(const EnergyPlusData &state,
+                  Real64 Eps, // required absolute accuracy
+                  int maxIters,
+                  int &SolFla,
+                  const std::function<Real64(Real64)> &f,
+                  Real64 X_0, // 1st bound of interval that contains the solution
+                  Real64 X_1, // 2nd bound of interval that contains the solution
+                  SolveRootStats &stats)
+{
+    // SUBROUTINE INFORMATION:
+    //       AUTHOR         Amir Roth
+    //       DATE WRITTEN   Nov. 2025
+
+    // PURPOSE OF THIS SUBROUTINE:
+    // This is a wrapper to SolveRoot that iterates over all root finding algorithms to find the best one.
+
+    Real64 XRes;
+
+    // Save and restore "global" root finding algorithm
+    RootAlgo algoTemp = state.dataRootFinder->rootAlgo;
+    state.dataRootFinder->rootAlgo = stats.algo;
+
+    SolveRoot(state, Eps, maxIters, SolFla, XRes, f, X_0, X_1);
+
+    state.dataRootFinder->rootAlgo = algoTemp;
+
+    if (SolFla > 0) {
+        stats.counts++;
+        stats.algoCounts[(int)stats.algo]++;
+        stats.algoIters[(int)stats.algo] += SolFla;
+
+        constexpr int TRIALS_PER_COUNT = 5;
+
+        // Trial period, cycle thru algorithms
+        if (stats.counts < TRIALS_PER_COUNT * (int)RootAlgo::Num) {
+            stats.algo = static_cast<RootAlgo>((int)stats.algo + 1);
+            if (stats.algo == RootAlgo::Num) {
+                stats.algo = RootAlgo::RegulaFalsi;
+            }
+
+            // Choose base algorithm, i.e., fewest total iterations
+        } else if (stats.counts == TRIALS_PER_COUNT * (int)RootAlgo::Num) {
+            int minIters = maxIters * TRIALS_PER_COUNT;
+            stats.algo = RootAlgo::Invalid;
+            for (int i = 0; i < (int)RootAlgo::Num; ++i) {
+                if (stats.algoIters[i] < minIters) {
+                    stats.algo = static_cast<RootAlgo>(i);
+                    minIters = stats.algoIters[i];
+                }
+            }
+
+            // Have chosen an algorithm, stats.algo should be it
+        } else {
+        }
+    }
+
+    return XRes;
+}
+
 void MovingAvg(Array1D<Real64> &DataIn, int const NumItemsInAvg)
 {
     if (NumItemsInAvg <= 1) {

src/EnergyPlus/General.hh

@@ -60,6 +60,8 @@
 #include <EnergyPlus/EPVector.hh>
 #include <EnergyPlus/EnergyPlus.hh>
 
+#include <EnergyPlus/HVACSystemRootFindingAlgorithm.hh>
+
 namespace EnergyPlus {
 
 // Forward declarations
@@ -72,6 +74,17 @@ namespace Weather {
 
 namespace General {
 
+    constexpr int SOLVEROOT_ERROR_INIT = -2;
+    constexpr int SOLVEROOT_ERROR_ITER = -1;
+
+    struct SolveRootStats
+    {
+        RootAlgo algo = RootAlgo::RegulaFalsi;
+        int counts = 0;
+        std::array<int, (int)RootAlgo::Num> algoCounts = {0};
+        std::array<int, (int)RootAlgo::Num> algoIters = {0};
+    };
+
     // A second version that does not require a payload -- use lambdas
     void SolveRoot(const EnergyPlusData &state,
                    Real64 Eps,   // required absolute accuracy
@@ -82,6 +95,15 @@ namespace General {
                    Real64 X_0,  // 1st bound of interval that contains the solution
                    Real64 X_1); // 2nd bound of interval that contains the solution
 
+    Real64 SolveRoot2(const EnergyPlusData &state,
+                      Real64 Eps,   // required absolute accuracy
+                      int maxIters, // maximum number of iterations
+                      int &SolFlag, // solution flag
+                      const std::function<Real64(Real64)> &f,
+                      Real64 X_0, // 1st bound of interval that contains the solution
+                      Real64 X_1,
+                      SolveRootStats &config); // 2nd bound of interval that contains the solution
+
     void MovingAvg(Array1D<Real64> &DataIn, int NumItemsInAvg);
 
     void ProcessDateString(EnergyPlusData &state,

src/EnergyPlus/HVACSystemRootFindingAlgorithm.hh

@@ -51,7 +51,8 @@
 #include <EnergyPlus/Data/BaseData.hh>
 
 namespace EnergyPlus {
-enum class HVACSystemRootSolverAlgorithm : int
+
+enum class RootAlgo : int
 {
     Invalid = -1,
     RegulaFalsi,
@@ -63,18 +64,16 @@ enum class HVACSystemRootSolverAlgorithm : int
     Num
 };
 
-static constexpr std::array<std::string_view, static_cast<int>(HVACSystemRootSolverAlgorithm::Num)> HVACSystemRootSolverAlgorithmUC = {
+static constexpr std::array<std::string_view, (int)RootAlgo::Num> rootAlgoNamesUC = {
     "REGULAFALSI", "BISECTION", "REGULAFALSITHENBISECTION", "BISECTIONTHENREGULAFALSI", "ALTERNATION", "SHORTBISECTIONTHENREGULAFALSI"};
 
-struct HVACSystemRootFindingAlgorithm
-{
-    std::string Algorithm = {}; // Choice of algorithm
-    int NumOfIter = 5;          // Number of Iteration Before Algorithm Switch
-    HVACSystemRootSolverAlgorithm HVACSystemRootSolverMethod = HVACSystemRootSolverAlgorithm::RegulaFalsi;
-};
 struct RootFindingData : BaseGlobalStruct
 {
-    HVACSystemRootFindingAlgorithm HVACSystemRootFinding;
+    std::string Algorithm = {}; // Choice of algorithm
+
+    int NumOfIter = 5; // Number of Iteration Before Algorith Switch
+    RootAlgo rootAlgo = RootAlgo::RegulaFalsi;
+
     void init_constant_state([[maybe_unused]] EnergyPlusData &state) override
     {
     }
@@ -85,7 +84,6 @@ struct RootFindingData : BaseGlobalStruct
 
     void clear_state() override
     {
-        this->HVACSystemRootFinding = {};
     }
 };
 

src/EnergyPlus/HeatBalanceManager.cc

@@ -513,9 +513,6 @@ namespace HeatBalanceManager {
         // ZoneAirHeatBalanceAlgorithm, Added by L. Gu, 12/09
         // ZoneAirContaminantBalance, Added by L. Gu, 06/10
 
-        // Using/Aliasing
-        auto &HVACSystemRootFinding = state.dataRootFinder->HVACSystemRootFinding;
-
         // SUBROUTINE PARAMETER DEFINITIONS:
         static constexpr std::string_view RoutineName("GetProjectControlData: ");
         static constexpr std::string_view routineName = "GetProjectControlData";
@@ -1224,31 +1221,30 @@ namespace HeatBalanceManager {
                                                                      state.dataIPShortCut->cNumericFieldNames);
             ErrorObjectHeader eoh{routineName, state.dataHeatBalMgr->CurrentModuleObject, AlphaName(1)};
             if (NumAlpha > 0) {
-                HVACSystemRootFinding.Algorithm = AlphaName(1);
-                HVACSystemRootFinding.HVACSystemRootSolverMethod =
-                    static_cast<HVACSystemRootSolverAlgorithm>(getEnumValue(HVACSystemRootSolverAlgorithmUC, Util::makeUPPER(AlphaName(1))));
-                if (HVACSystemRootFinding.HVACSystemRootSolverMethod == HVACSystemRootSolverAlgorithm::Invalid) {
-                    HVACSystemRootFinding.HVACSystemRootSolverMethod = HVACSystemRootSolverAlgorithm::RegulaFalsi;
+                state.dataRootFinder->Algorithm = AlphaName(1);
+                state.dataRootFinder->rootAlgo = static_cast<RootAlgo>(getEnumValue(rootAlgoNamesUC, Util::makeUPPER(AlphaName(1))));
+                if (state.dataRootFinder->rootAlgo == RootAlgo::Invalid) {
+                    state.dataRootFinder->rootAlgo = RootAlgo::RegulaFalsi;
                     ShowWarningInvalidKey(
                         state, eoh, state.dataIPShortCut->cAlphaFieldNames(1), AlphaName(1), "Invalid input. The default choice is assigned.");
                     ShowContinueError(
                         state, "Valid choices are: RegulaFalsi, Bisection, BisectionThenRegulaFalsi, RegulaFalsiThenBisection, or Alternation.");
                 }
             }
             if (NumNumber > 0) {
-                HVACSystemRootFinding.NumOfIter = BuildingNumbers(1);
+                state.dataRootFinder->NumOfIter = BuildingNumbers(1);
             }
         } else {
-            HVACSystemRootFinding.Algorithm = "RegulaFalsi";
-            HVACSystemRootFinding.HVACSystemRootSolverMethod = HVACSystemRootSolverAlgorithm::RegulaFalsi;
+            state.dataRootFinder->Algorithm = "RegulaFalsi";
+            state.dataRootFinder->rootAlgo = RootAlgo::RegulaFalsi;
         }
 
         // Write Solution Algorithm to the initialization output file for User Verification
         constexpr const char *Format_734(
             "! <HVACSystemRootFindingAlgorithm>, Value {{RegulaFalsi | Bisection | BisectionThenRegulaFalsi | RegulaFalsiThenBisection}}\n");
         constexpr const char *Format_735(" HVACSystemRootFindingAlgorithm, {}\n");
         print(state.files.eio, Format_734);
-        print(state.files.eio, Format_735, HVACSystemRootFinding.Algorithm);
+        print(state.files.eio, Format_735, state.dataRootFinder->Algorithm);
     }
 
     void GetSiteAtmosphereData(EnergyPlusData &state, bool &ErrorsFound)