Fix #11867 - remove boost::bisect and correct calculation of Standard Effective Temperature by jmarrec · Pull Request #11877 · NatLabRockies/EnergyPlus

jmarrec · 2026-10-23T15:20:45Z

Pull request overview

Fixes boost::math::tools::bisect is used in E+ own code #11867

REVIEW WITH https://github.com/NatLabRockies/EnergyPlus/pull/11877/changes/29ec00e6e1e616328f7fc6eb554b21d2db6f1dac..7075e054ccb57652ac8d2eb635e9f406bb86fc6f . The first commit just takes an alias for readability.

Description of the purpose of this PR

Initially I changed boost::bisect to SolveRoot and found that SET is discontinuous.
Fix SET discontinuity: correct SkinWetSweat when skin wettedness is capped at EvapEff
- We now deliberately differ from the ASHRAE 55 reference code

TL;DR is that the reference implementation and tables from ASRHAE 55 Appendices have an error (that the original Gagge et al. (1986) didn't have). Wasn't sure how to approach this, so I reported it to, and then fixed myself two open source implementations that are probably much closer to ASHRAE 55 (I suspect at least of the two is actually the reference...):

Should watch:

ASHRAE 55: submit an addendum correcting p_rsw at the wettedness cap in the SET code pythermalcomfort/pythermalcomfort#492

Pull Request Author

Reviewer

…adable...

…ffect SET is discontinuous at the critical skin wettedness cap, so SolveRoot, which converges on |f(x)|, failed on 42 timesteps of 1ZoneUncontrolled_Win_ASH55_Thermal_Comfort. Bisecting on the bracket width matches boost::math::tools::bisect, and the results are identical to develop. cf: - #11867 (comment)

…mperature (SET) ``` [ RUN ] EnergyPlusFixture.ThermalComfort_CalcStandardEffectiveTemp_ContinuousAtCriticalWettedness /media/DataExt4/Software/Others/EnergyPlus/tst/EnergyPlus/unit/ThermalComfort.unit.cc:1133: Failure Expected: (maxStep) < (0.01), actual: 0.486617 vs 0.01 /media/DataExt4/Software/Others/EnergyPlus/tst/EnergyPlus/unit/ThermalComfort.unit.cc:1134: Failure Expected: (minStep) > (0.0), actual: -0.486617 vs 0 [ FAILED ] EnergyPlusFixture.ThermalComfort_CalcStandardEffectiveTemp_ContinuousAtCriticalWettedness (607 ms) ```

…ness is capped at EvapEff SkinWetSweat = EvapEff / 0.94 is not the inverse of SkinWetTot = 0.06 + 0.94 * SkinWetSweat, so at the cap the skin evaporated as if its wettedness were EvapEff + 0.06. The cap then switched on and off between the minutes of the simulation, making SET jump by up to 0.5 C. Use (SkinWetTot - 0.06) / 0.94, as in Gagge et al. (1986): prsw=(eveff-.06)/.94. Same fix as: - CenterForTheBuiltEnvironment/comfort_tool#103 - pythermalcomfort/pythermalcomfort#490

…continuous

SolveRoot bounds the SET mismatch, not the cooling effect itself, so Acc = 0.01 allowed errors up to 0.017 C where SET rises slowly with temperature (vs 0.003 C for the previous bisection). With 0.001, the max error on the annual ASH55 inputs is 0.0017 C for 4.4 root function evaluations per call on average, vs 15 with bisection.

jmarrec · 2026-10-23T15:25:02Z

            if (s_tc->SkinWetTot >= EvapEff) {
                s_tc->SkinWetTot = EvapEff;
-                s_tc->SkinWetSweat = EvapEff / 0.94;
+                // Invert SkinWetTot = SkinWetSweat + SkinWetDiff = SkinWetSweat + (1 - SkinWetSweat) * 0.06 = 0.06 + 0.94 * SkinWetSweat, so
+                // that EvapHeatLoss below equals SkinWetTot * EvapHeatLossMax (Gagge et al. 1986: prsw=(eveff-.06)/.94)
+                s_tc->SkinWetSweat = (s_tc->SkinWetTot - 0.06) / 0.94;


That's the discontinuity fix.

jmarrec · 2026-10-23T15:28:56Z

+        Real64 constexpr tempAccuracy = 0.001; // Accuracy of the SET match [deltaC]
+        int constexpr MaxIte = 100;
+        Real64 constexpr lowerBound = 0.0;
+        Real64 constexpr upperBound = 50.0;
+        int SolFla = 0;

-        // We have yet another solver?
-        try {
-            std::pair<Real64, Real64> solverResult = boost::math::tools::bisect(ce_root_function, lowerBound, upperBound, ce_root_termination);
-            CoolingEffect = (solverResult.first + solverResult.second) / 2;
-        } catch (const std::exception &e) {
+        General::SolveRoot(state, tempAccuracy, MaxIte, SolFla, CoolingEffect, ce_root_function, lowerBound, upperBound);
+        if (SolFla < 0) {


Now that SET is continuous, we can go back to General::SolveRoot instead of the bisection.

I first used tempAccuracy = 0.01 °C, but then lowered it to 0.001. SolveRoot stops when the SET mismatch |f(x)| is below the tolerance, not when the bracket on the cooling effect is small. Where SET rises slowly with temperature, a 0.01 tolerance on SET can therefore allow a larger error on the cooling effect.

I instrumented a standalone copy of the solver and ran it on the 38016 (air temperature, MRT, RH, air velocity, met, clo) inputs from the annual ASH55 run. I compared each result with the exact root:

Solver Max error [°C] Mean error [°C] Root function evaluations per call

Bisection (previous) 0.0031 0.0015 15 (always)

SolveRoot, Acc = 0.01 0.0169 0.0012 3.98

SolveRoot, Acc = 0.001 0.0017 0.0001 4.39

With 0.001, SolveRoot is more accurate than the bisection. It still needs about 3.4x fewer evaluations, and each evaluation is a full SET calculation. The bisection always does 15, because it stops on bracket width: two endpoints, then halvings until the [0, 50] bracket is under 0.01.

The annual E+ run (fix + bisection against fix + SolveRoot) shows the same thing. The max Cooling Effect diff is 0.018 °C with Acc = 0.01 and 0.004 °C with Acc = 0.001.

Improved accuracy with fewer iterations = winning

(I believe the note above where you state "3.4x" should probably be 4.3x based on the table.)

jmarrec · 2026-10-23T15:29:24Z

+TEST_F(EnergyPlusFixture, ThermalComfort_CalcStandardEffectiveTemp_ContinuousAtCriticalWettedness)
+{
+    // When the skin wettedness reaches EvapEff, SkinWetSweat must invert SkinWetTot = 0.06 + 0.94 * SkinWetSweat, i.e.
+    // (SkinWetTot - 0.06) / 0.94 as in Gagge et al. (1986). SET should then change smoothly and monotonically with temperature.
+    state->dataEnvrn->OutBaroPress = 101325.0;
+    Real64 constexpr RelHum = 0.7;
+    Real64 constexpr AirVel = 0.1;
+    Real64 constexpr ActMet = 1.1;
+    Real64 constexpr CloUnit = 0.5;
+    Real64 constexpr WorkEff = 0.0;
+
+    Real64 maxStep = 0.0;
+    Real64 minStep = std::numeric_limits<Real64>::max();
+    Real64 previous = CalcStandardEffectiveTemp(*state, 32.5, 32.5, RelHum, AirVel, ActMet, CloUnit, WorkEff);
+    for (int i = 1; i <= 1500; ++i) {
+        Real64 const temp = 32.5 + i * 0.001;
+        Real64 const current = CalcStandardEffectiveTemp(*state, temp, temp, RelHum, AirVel, ActMet, CloUnit, WorkEff);
+        maxStep = std::max(maxStep, std::abs(current - previous));
+        minStep = std::min(minStep, current - previous);
+        previous = current;
+    }
+    EXPECT_LT(maxStep, 0.01);
+    EXPECT_GT(minStep, 0.0); // monotonic, so the cooling effect has a single root
+}


Before fix

[ RUN ] EnergyPlusFixture.ThermalComfort_CalcStandardEffectiveTemp_ContinuousAtCriticalWettedness /media/DataExt4/Software/Others/EnergyPlus/tst/EnergyPlus/unit/ThermalComfort.unit.cc:1133: Failure Expected: (maxStep) < (0.01), actual: 0.486617 vs 0.01 /media/DataExt4/Software/Others/EnergyPlus/tst/EnergyPlus/unit/ThermalComfort.unit.cc:1134: Failure Expected: (minStep) > (0.0), actual: -0.486617 vs 0 [ FAILED ] EnergyPlusFixture.ThermalComfort_CalcStandardEffectiveTemp_ContinuousAtCriticalWettedness (607 ms)

Looks good. There's a few lingering test values in here where @Myoldmopar was scratching his head previously.

https://github.com/NatLabRockies/EnergyPlus/blame/7075e054ccb57652ac8d2eb635e9f406bb86fc6f/tst/EnergyPlus/unit/ThermalComfort.unit.cc#L1085

Might be worth cleaning these up while we're in here.

github-actions · 2026-10-23T15:38:44Z

⚠️ Regressions detected on macos-15 for commit `62778fa`

Regression Summary

ESO Big Diffs: 4

jmarrec · 2026-10-23T15:45:02Z

Diffs

Expected ESO diffs in 4 files, and no table diffs:

File	Max SET diff [°C]	Max Cooling Effect diff [°C]
`1ZoneUncontrolled_Win_ASH55_Thermal_Comfort`	+0.270	0.006
`1ZoneUncontrolled_Win_ResilienceReports`	+0.283	—
`1ZoneUncontrolled_Win_ResilienceReports_ReportingPeriod`	+0.281	—
`1ZoneUncontrolled_Win_ResilienceReports_ReportingPeriod_AllTable`	+0.281	—

These are the only files that report SET (and, for ASH55, the Cooling Effect). All the max diffs are at 07/21 06:00 on the summer design day. The Cooling Effect Adjusted PMV/PPD diffs in ASH55 are tiny (max 0.0023 and 0.013).

Annual run

To check for a pattern, I ran 1ZoneUncontrolled_Win_ASH55_Thermal_Comfort for a full year (Golden, CO TMY3) at timestep resolution (35040 timesteps), on develop and on this branch.
I changed the variables from Hourly to Timestep.

Only the SET, Cooling Effect and Cooling Effect Adjusted PMV/PPD outputs change.

SET changes at 7550 timesteps, all with SET above 36.6 °C, which is where skin wettedness reaches the cap. Below that, the two builds are bit-identical. The change is a systematic upward correction. With the bug, the skin evaporated an extra 0.06 * EvapHeatLossMax, so SET came out too low.

SET (develop) [°C]	Mean change [°C]	Max change [°C]
37–39	+0.31	+0.49
41–42	+0.24	+0.34
43–44	+0.16	+0.24
45–46	+0.09	+0.15
47–48	+0.02	+0.05

The correction shrinks as conditions get hotter and more humid, because EvapHeatLossMax, and with it the extra 0.06 * EvapHeatLossMax, gets smaller. There are 14 timesteps around 37–38 °C where SET decreases slightly, by at most 0.003 °C. At those timesteps the cap was toggling on and off.

Cooling Effect:

At timesteps where SET changed, the Cooling Effect follows. At 678 timesteps (1.9 % of the year, about 4.6 % from June to September), it changes by 0.15 to 0.47 °C, nearly always downward (672 of 678). On develop, the SET discontinuity gave f(x) = SET(still air, T - x) - SET(v, T) several roots, and the solver landed on a wrong one. SET is now continuous and monotonic, so there is a single root.
At timesteps where SET is unchanged, the only differences come from switching the bisection to SolveRoot: max 0.0076 °C, mean 0.0015 °C, p99 0.0032 °C.

The largest Cooling Effect Adjusted PMV and PPD changes are 0.13 and 4.7, at the wrong-root timesteps.

No "The cooling effect could not be solved" warnings in either run.

github-actions · 2026-10-23T15:46:28Z

⚠️ Regressions detected on ubuntu-24.04 for commit `62778fa`

Regression Summary

ESO Big Diffs: 4

mitchute

Thanks for all the background work to figure this out @jmarrec. Looks good. I'll let you decide if you want to address those minor cleanups in the unit tests.

mitchute · 2026-10-29T04:30:21Z

-// Third Party Headers
-#include <boost/math/tools/roots.hpp>


That's great.

mitchute · 2026-10-29T04:30:50Z

    Real64 CalcStandardEffectiveTemp(
        EnergyPlusData &state, Real64 AirTemp, Real64 RadTemp, Real64 RelHum, Real64 AirVel, Real64 ActMet, Real64 CloUnit, Real64 WorkEff)
    {
+        auto &s_tc = state.dataThermalComforts;


More comfortable to read.

mitchute · 2026-10-29T04:31:49Z

            if (s_tc->SkinWetTot >= EvapEff) {
                s_tc->SkinWetTot = EvapEff;
-                s_tc->SkinWetSweat = EvapEff / 0.94;
+                // Invert SkinWetTot = SkinWetSweat + SkinWetDiff = SkinWetSweat + (1 - SkinWetSweat) * 0.06 = 0.06 + 0.94 * SkinWetSweat, so
+                // that EvapHeatLoss below equals SkinWetTot * EvapHeatLossMax (Gagge et al. 1986: prsw=(eveff-.06)/.94)
+                s_tc->SkinWetSweat = (s_tc->SkinWetTot - 0.06) / 0.94;


mitchute · 2026-10-29T04:35:08Z

+        Real64 constexpr tempAccuracy = 0.001; // Accuracy of the SET match [deltaC]
+        int constexpr MaxIte = 100;
+        Real64 constexpr lowerBound = 0.0;
+        Real64 constexpr upperBound = 50.0;
+        int SolFla = 0;

-        // We have yet another solver?
-        try {
-            std::pair<Real64, Real64> solverResult = boost::math::tools::bisect(ce_root_function, lowerBound, upperBound, ce_root_termination);
-            CoolingEffect = (solverResult.first + solverResult.second) / 2;
-        } catch (const std::exception &e) {
+        General::SolveRoot(state, tempAccuracy, MaxIte, SolFla, CoolingEffect, ce_root_function, lowerBound, upperBound);
+        if (SolFla < 0) {


Improved accuracy with fewer iterations = winning

(I believe the note above where you state "3.4x" should probably be 4.3x based on the table.)

mitchute · 2026-10-29T04:38:32Z

+TEST_F(EnergyPlusFixture, ThermalComfort_CalcStandardEffectiveTemp_ContinuousAtCriticalWettedness)
+{
+    // When the skin wettedness reaches EvapEff, SkinWetSweat must invert SkinWetTot = 0.06 + 0.94 * SkinWetSweat, i.e.
+    // (SkinWetTot - 0.06) / 0.94 as in Gagge et al. (1986). SET should then change smoothly and monotonically with temperature.
+    state->dataEnvrn->OutBaroPress = 101325.0;
+    Real64 constexpr RelHum = 0.7;
+    Real64 constexpr AirVel = 0.1;
+    Real64 constexpr ActMet = 1.1;
+    Real64 constexpr CloUnit = 0.5;
+    Real64 constexpr WorkEff = 0.0;
+
+    Real64 maxStep = 0.0;
+    Real64 minStep = std::numeric_limits<Real64>::max();
+    Real64 previous = CalcStandardEffectiveTemp(*state, 32.5, 32.5, RelHum, AirVel, ActMet, CloUnit, WorkEff);
+    for (int i = 1; i <= 1500; ++i) {
+        Real64 const temp = 32.5 + i * 0.001;
+        Real64 const current = CalcStandardEffectiveTemp(*state, temp, temp, RelHum, AirVel, ActMet, CloUnit, WorkEff);
+        maxStep = std::max(maxStep, std::abs(current - previous));
+        minStep = std::min(minStep, current - previous);
+        previous = current;
+    }
+    EXPECT_LT(maxStep, 0.01);
+    EXPECT_GT(minStep, 0.0); // monotonic, so the cooling effect has a single root
+}


Looks good. There's a few lingering test values in here where @Myoldmopar was scratching his head previously.

https://github.com/NatLabRockies/EnergyPlus/blame/7075e054ccb57652ac8d2eb635e9f406bb86fc6f/tst/EnergyPlus/unit/ThermalComfort.unit.cc#L1085

Might be worth cleaning these up while we're in here.

jmarrec added 7 commits October 22, 2026 17:35

Take an alias in CalcStandardEffectiveTemp just so it's a bit more re…

29ec00e

…adable...

Fix #11867 - Replace boost bisect with SolveRoot

fad44d2

#11867 - Go back to SolveRoot for the cooling effect now that SET is …

1628097

…continuous

jmarrec self-assigned this Oct 23, 2026

jmarrec added the Defect Includes code to repair a defect in EnergyPlus label Oct 23, 2026

jmarrec commented Oct 23, 2026

View reviewed changes

jmarrec requested a review from mitchute October 26, 2026 18:42

mitchute approved these changes Oct 29, 2026

View reviewed changes

mitchute added the Staged By adding the "Staged" label, some expensive checks needed before merge are run. label Oct 29, 2026

Solver	Max error [°C]	Mean error [°C]	Root function evaluations per call
Bisection (previous)	0.0031	0.0015	15 (always)
SolveRoot, Acc = 0.01	0.0169	0.0012	3.98
SolveRoot, Acc = 0.001	0.0017	0.0001	4.39

Conversation

jmarrec commented Oct 23, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Pull request overview

Description of the purpose of this PR

Pull Request Author

Reviewer

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

github-actions Bot commented Oct 23, 2026

⚠️ Regressions detected on macos-15 for commit 62778fa

Uh oh!

jmarrec commented Oct 23, 2026

Diffs

Annual run

Uh oh!

github-actions Bot commented Oct 23, 2026

⚠️ Regressions detected on ubuntu-24.04 for commit 62778fa

Uh oh!

mitchute left a comment

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

mitchute Oct 29, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

This branch has not been deployed

Reviewers

Uh oh!

Assignees

Labels

Projects

Milestone

Development

Uh oh!

3 participants

jmarrec commented Oct 23, 2026 •

edited

Loading

⚠️ Regressions detected on macos-15 for commit `62778fa`

⚠️ Regressions detected on ubuntu-24.04 for commit `62778fa`

mitchute Oct 29, 2026 •

edited

Loading