Line data Source code
1 : // EnergyPlus, Copyright (c) 1996-2024, The Board of Trustees of the University of Illinois,
2 : // The Regents of the University of California, through Lawrence Berkeley National Laboratory
3 : // (subject to receipt of any required approvals from the U.S. Dept. of Energy), Oak Ridge
4 : // National Laboratory, managed by UT-Battelle, Alliance for Sustainable Energy, LLC, and other
5 : // contributors. All rights reserved.
6 : //
7 : // NOTICE: This Software was developed under funding from the U.S. Department of Energy and the
8 : // U.S. Government consequently retains certain rights. As such, the U.S. Government has been
9 : // granted for itself and others acting on its behalf a paid-up, nonexclusive, irrevocable,
10 : // worldwide license in the Software to reproduce, distribute copies to the public, prepare
11 : // derivative works, and perform publicly and display publicly, and to permit others to do so.
12 : //
13 : // Redistribution and use in source and binary forms, with or without modification, are permitted
14 : // provided that the following conditions are met:
15 : //
16 : // (1) Redistributions of source code must retain the above copyright notice, this list of
17 : // conditions and the following disclaimer.
18 : //
19 : // (2) Redistributions in binary form must reproduce the above copyright notice, this list of
20 : // conditions and the following disclaimer in the documentation and/or other materials
21 : // provided with the distribution.
22 : //
23 : // (3) Neither the name of the University of California, Lawrence Berkeley National Laboratory,
24 : // the University of Illinois, U.S. Dept. of Energy nor the names of its contributors may be
25 : // used to endorse or promote products derived from this software without specific prior
26 : // written permission.
27 : //
28 : // (4) Use of EnergyPlus(TM) Name. If Licensee (i) distributes the software in stand-alone form
29 : // without changes from the version obtained under this License, or (ii) Licensee makes a
30 : // reference solely to the software portion of its product, Licensee must refer to the
31 : // software as "EnergyPlus version X" software, where "X" is the version number Licensee
32 : // obtained under this License and may not use a different name for the software. Except as
33 : // specifically required in this Section (4), Licensee shall not use in a company name, a
34 : // product name, in advertising, publicity, or other promotional activities any name, trade
35 : // name, trademark, logo, or other designation of "EnergyPlus", "E+", "e+" or confusingly
36 : // similar designation, without the U.S. Department of Energy's prior written consent.
37 : //
38 : // THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR
39 : // IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY
40 : // AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR
41 : // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
42 : // CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
43 : // SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
44 : // THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR
45 : // OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
46 : // POSSIBILITY OF SUCH DAMAGE.
47 :
48 : // C++ Headers
49 : #include <cmath>
50 :
51 : // ObjexxFCL Headers
52 : #include <ObjexxFCL/Fmath.hh>
53 :
54 : // EnergyPlus Headers
55 : #include <EnergyPlus/Data/EnergyPlusData.hh>
56 : #include <EnergyPlus/DataPrecisionGlobals.hh>
57 : #include <EnergyPlus/General.hh>
58 : #include <EnergyPlus/RootFinder.hh>
59 : #include <EnergyPlus/UtilityRoutines.hh>
60 :
61 : namespace EnergyPlus::RootFinder {
62 : // Module containing the derivative-free, root finder routines for smooth,
63 : // "mostly" monotonic functions in one variable.
64 : // The algorithm is guaranteed to find the root if:
65 : // - it is bracketed between min/max bounds
66 : // - there is an odd number of roots between the min/max bounds
67 : // Note that there is an even number of roots between the min/max bounds then it is possible
68 : // that the algorithm will terminate with a min or max constrained solution instead of the
69 : // actual root.
70 :
71 : // MODULE INFORMATION:
72 : // AUTHOR Dimitri Curtil (LBNL)
73 : // DATE WRITTEN February 2006
74 : // MODIFIED na
75 : // RE-ENGINEERED na
76 :
77 : // PURPOSE OF THIS MODULE:
78 : // To encapsulate the data and algorithms required to
79 : // find the root of a 1-dimensional function F(X) = Y .
80 : // At any moment during the solution process, the following condition is satisfied:
81 : // - For an overall increasing function F(X):
82 : // MinPoint%X <= LowerPoint%X <= CurrentPoint%X <= UpperPoint%X <= MaxPoint%X
83 : // MinPoint%Y <= LowerPoint%Y <= CurrentPoint%Y <= UpperPoint%Y <= MaxPoint%Y
84 : // - For an overall decreasing function F(X):
85 : // MinPoint%X <= LowerPoint%X <= CurrentPoint%X <= UpperPoint%X <= MaxPoint%X
86 : // MinPoint%Y >= LowerPoint%Y >= CurrentPoint%Y >= UpperPoint%Y >= MaxPoint%Y
87 : // Requirements for function F():
88 : // - smooth, continuous, deterministic function
89 : // - root XRoot bracketed within the min and max bounds, referred to as XMin, XMax
90 : // - XMin and XMax must remain constant during iterations
91 : // - If F(X) is increasing, set action to iSlopeIncreasing
92 : // This implies that if X1 < X2, then F(X1) < F(X2) must be satisfied.
93 : // - If F(X) is decreasing, set action to iSlopeDecreasing
94 : // This implies that if X1 < X2, then F(X1) > F(X2) must be satisfied.
95 : // Note that the residual function is allowed to be non-monotonic between the min and max
96 : // points as long as the slope specification is satisfied at the min and max points.
97 : // Unconstrained convergence conditions:
98 : // - ABS(F(XRoot)) <= ATolY
99 : // - ABS(X2-X1) <= TolX * (ABS(X2)+ABS(X1))/2 + ATolX
100 : // Constrained convergence conditions:
101 : // For iSlopeIncreasing action:
102 : // - YMin >= 0
103 : // - YMax <= 0
104 : // For iSlopeDecreasing action:
105 : // - YMin <= 0
106 : // - YMax >= 0
107 :
108 : // METHODOLOGY EMPLOYED:
109 : // Implements an hybrid solution method consisting of:
110 : // - bisection method (aka interval halving)
111 : // - false position method (aka regula falsi)
112 : // - secant method
113 : // - Brent's method (inverse quadratic interpolation)
114 : // with safeguards against:
115 : // - out-of-range error
116 : // - singularity (i.e., locally singular Jacobian)
117 : // - bad slope
118 : // - round-off error
119 : // Typical usage of the root finder, assuming that :
120 : // - the function is called MyFunc()
121 : // - the function is strictly increasing between the min and max bounds
122 : // - the root is bracketed between XMin=0 and XMax=1
123 : // - the function is defined with the prototype real FUNCTION MyFunc( X )
124 : // - the solution method to use is the Brent's method
125 : // - the absolute tolerance for Y=MyFunc(X) is ATolY=1.0E-4
126 : // As a safeguard the iterative process stops after 50 iterations if the root has not been
127 : // located within the specified tolerance yet.
128 : // <BEGIN CODE SNIPPET>
129 : // TYPE(RootFinderDataType) :: RF
130 : // INTEGER :: IterationCount
131 : // REAL(r64) :: X, Y
132 : // LOGICAL :: IsDoneFlag
133 : // CALL SetupRootFinder( &
134 : // RF, & ! RootFinderData
135 : // iSlopeIncreasing, & ! SlopeType
136 : // iMethodBrent, & ! MethodType
137 : // 1.0d-6, & ! TolX
138 : // 1.0d-6, & ! ATolX
139 : // 1.0d-4 & ! ATolY
140 : // )
141 : // CALL InitializeRootFinder( &
142 : // RF, & ! RootFinderData
143 : // 0.0d0, & ! XMin
144 : // 1.0d0 & ! XMax
145 : // )
146 : // IterationCount = 0
147 : // IsDoneFlag = .FALSE.
148 : // RF%XCandidate = 0.1d0 ! Sets X to initial value XInit
149 : // DO WHILE ( .NOT.IsDoneFlag )
150 : // IterationCount = IterationCount+1
151 : // IF ( IterationCount>50 ) THEN
152 : // WRITE(*,*) 'Error: Too many iterations..."
153 : // EXIT
154 : // END IF
155 : // ! Evaluate function with new root candidate
156 : // X = RF%XCandidate
157 : // Y = MyFunc( X )
158 : // ! Update root finder data with new iterate
159 : // CALL IterateRootFinder( RF, X, Y, IsDoneFlag )
160 : // END DO
161 : // ! Write root finder status description to standard output
162 : // CALL WriteRootFinderStatus( 6, RF )
163 : // <END CODE SNIPPET>
164 :
165 : // REFERENCES:
166 : // "Numerical Recipes in FORTRAN", Chapter 9 "Root Finding and Nonlinear Sets of Equations", pp.340-352
167 : // Used for formulas, but not for code.
168 :
169 : // Using/Aliasing
170 : using namespace DataRootFinder;
171 :
172 125006 : void SetupRootFinder(EnergyPlusData &state,
173 : RootFinderDataType &RootFinderData, // Data used by root finding algorithm
174 : DataRootFinder::Slope const SlopeType, // Either Slope::Increasing or Slope::Decreasing
175 : DataRootFinder::RootFinderMethod const MethodType, // Any of the iMethod<name> code but iMethodNone
176 : Real64 const TolX, // Relative tolerance for X variables
177 : Real64 const ATolX, // Absolute tolerance for X variables
178 : Real64 const ATolY // Absolute tolerance for Y variables
179 : )
180 : {
181 :
182 : // SUBROUTINE INFORMATION:
183 : // AUTHOR Dimitri Curtil (LBNL)
184 : // DATE WRITTEN February 2006
185 : // MODIFIED
186 : // RE-ENGINEERED na
187 :
188 : // PURPOSE OF THIS SUBROUTINE:
189 : // This subroutine loads the numerical controls for the root finder.
190 :
191 : // Load assumed action for underlying function F(X)
192 125006 : if (SlopeType != DataRootFinder::Slope::Increasing && SlopeType != DataRootFinder::Slope::Decreasing) {
193 0 : ShowSevereError(state, "SetupRootFinder: Invalid function slope specification. Valid choices are:");
194 0 : ShowContinueError(state, format("SetupRootFinder: Slope::Increasing={}", DataRootFinder::Slope::Increasing));
195 0 : ShowContinueError(state, format("SetupRootFinder: Slope::Decreasing={}", DataRootFinder::Slope::Decreasing));
196 0 : ShowFatalError(state, "SetupRootFinder: Preceding error causes program termination.");
197 : }
198 125006 : RootFinderData.Controls.SlopeType = SlopeType;
199 :
200 : // Load solution method
201 125006 : if (MethodType != RootFinderMethod::Bisection && MethodType != RootFinderMethod::FalsePosition && MethodType != RootFinderMethod::Secant &&
202 : MethodType != RootFinderMethod::Brent) {
203 :
204 0 : ShowSevereError(state, "SetupRootFinder: Invalid solution method specification. Valid choices are:");
205 0 : ShowContinueError(state, format("SetupRootFinder: iMethodBisection={}", RootFinderMethod::Bisection));
206 0 : ShowContinueError(state, format("SetupRootFinder: iMethodFalsePosition={}", RootFinderMethod::FalsePosition));
207 0 : ShowContinueError(state, format("SetupRootFinder: iMethodSecant={}", RootFinderMethod::Secant));
208 0 : ShowContinueError(state, format("SetupRootFinder: iMethodBrent={}", RootFinderMethod::Brent));
209 0 : ShowFatalError(state, "SetupRootFinder: Preceding error causes program termination.");
210 : }
211 125006 : RootFinderData.Controls.MethodType = MethodType;
212 :
213 : // Load relative tolerance parameter for X variables
214 125006 : if (TolX < 0.0) {
215 0 : ShowFatalError(state, "SetupRootFinder: Invalid tolerance specification for X variables. TolX >= 0");
216 : }
217 125006 : RootFinderData.Controls.TolX = TolX;
218 :
219 : // Load absolute tolerance parameter for X variables
220 125006 : if (ATolX < 0.0) {
221 0 : ShowFatalError(state, "SetupRootFinder: Invalid absolute tolerance specification for X variables. ATolX >= 0");
222 : }
223 125006 : RootFinderData.Controls.ATolX = ATolX;
224 :
225 : // Load absolute tolerance parameter for Y variables
226 125006 : if (ATolY < 0.0) {
227 0 : ShowFatalError(state, "SetupRootFinder: Invalid absolute tolerance specification for Y variables. ATolY >= 0");
228 : }
229 125006 : RootFinderData.Controls.ATolY = ATolY;
230 :
231 : // Reset internal data for root finder with fictive min and max values
232 125006 : ResetRootFinder(RootFinderData, DataPrecisionGlobals::constant_zero, DataPrecisionGlobals::constant_zero);
233 125006 : }
234 :
235 9047304 : void ResetRootFinder(RootFinderDataType &RootFinderData, // Data used by root finding algorithm
236 : Real64 const XMin, // Minimum X value allowed
237 : Real64 const XMax // Maximum X value allowed
238 : )
239 : {
240 :
241 : // SUBROUTINE INFORMATION:
242 : // AUTHOR Dimitri Curtil (LBNL)
243 : // DATE WRITTEN February 2006
244 : // MODIFIED
245 : // RE-ENGINEERED na
246 :
247 : // PURPOSE OF THIS SUBROUTINE:
248 : // This subroutine initializes the data for the root finder.
249 :
250 : // Reset min point
251 9047304 : RootFinderData.MinPoint.X = XMin;
252 9047304 : RootFinderData.MinPoint.Y = 0.0;
253 9047304 : RootFinderData.MinPoint.DefinedFlag = false;
254 :
255 : // Reset max point
256 9047304 : RootFinderData.MaxPoint.X = XMax;
257 9047304 : RootFinderData.MaxPoint.Y = 0.0;
258 9047304 : RootFinderData.MaxPoint.DefinedFlag = false;
259 :
260 : // Reset lower point
261 9047304 : RootFinderData.LowerPoint.X = 0.0;
262 9047304 : RootFinderData.LowerPoint.Y = 0.0;
263 9047304 : RootFinderData.LowerPoint.DefinedFlag = false;
264 :
265 : // Reset upper point
266 9047304 : RootFinderData.UpperPoint.X = 0.0;
267 9047304 : RootFinderData.UpperPoint.Y = 0.0;
268 9047304 : RootFinderData.UpperPoint.DefinedFlag = false;
269 :
270 : // Reset previous point
271 9047304 : RootFinderData.CurrentPoint.X = 0.0;
272 9047304 : RootFinderData.CurrentPoint.Y = 0.0;
273 9047304 : RootFinderData.CurrentPoint.DefinedFlag = false;
274 :
275 : // Reset iterate history with last 3 best points
276 9047304 : RootFinderData.NumHistory = 0;
277 36189216 : for (auto &e : RootFinderData.History) {
278 27141912 : e.X = e.Y = 0.0;
279 27141912 : e.DefinedFlag = false;
280 : }
281 :
282 : // Reset increments over successive iterations
283 9047304 : RootFinderData.Increment.X = 0.0;
284 9047304 : RootFinderData.Increment.Y = 0.0;
285 9047304 : RootFinderData.Increment.DefinedFlag = false;
286 :
287 9047304 : RootFinderData.XCandidate = 0.0;
288 :
289 : // Reset default state
290 9047304 : RootFinderData.StatusFlag = RootFinderStatus::None;
291 9047304 : RootFinderData.CurrentMethodType = RootFinderMethod::None;
292 9047304 : RootFinderData.ConvergenceRate = -1.0;
293 9047304 : }
294 :
295 8922298 : void InitializeRootFinder(EnergyPlusData &state,
296 : RootFinderDataType &RootFinderData, // Data used by root finding algorithm
297 : Real64 const XMin, // Minimum X value allowed
298 : Real64 const XMax // Maximum X value allowed
299 : )
300 : {
301 :
302 : // SUBROUTINE INFORMATION:
303 : // AUTHOR Dimitri Curtil (LBNL)
304 : // DATE WRITTEN February 2006
305 : // MODIFIED
306 : // RE-ENGINEERED na
307 :
308 : // PURPOSE OF THIS SUBROUTINE:
309 : // This subroutine initializes the min and max for the root finder before
310 : // finding a new root.
311 :
312 : // SUBROUTINE LOCAL VARIABLE DECLARATIONS:
313 : Real64 SavedXCandidate;
314 : Real64 XMinReset;
315 :
316 8922298 : XMinReset = XMin;
317 8922298 : if (XMin > XMax) {
318 0 : if (XMax == 0.0) {
319 0 : XMinReset = XMax;
320 : } else {
321 0 : ShowFatalError(state, format("InitializeRootFinder: Invalid min/max bounds XMin={:.6T} must be smaller than XMax={:.6T}", XMin, XMax));
322 : }
323 : }
324 :
325 : // First save current candidate value before it is overriden in ResetRootFinder()
326 8922298 : SavedXCandidate = RootFinderData.XCandidate;
327 :
328 : // Reset internal data for root finder with actual min and max values
329 : // NOTE: This resets the value of RootFinderData%XCandidate to zero
330 8922298 : ResetRootFinder(RootFinderData, XMinReset, XMax);
331 :
332 : // Enforce min/max constraints on previous candidate if available
333 : // NOTE: If XMin == XMax then this forces the value of XCandidateto the desired solution
334 8922298 : RootFinderData.XCandidate = min(RootFinderData.MaxPoint.X, max(SavedXCandidate, RootFinderData.MinPoint.X));
335 8922298 : }
336 :
337 19745373 : void IterateRootFinder(EnergyPlusData &state,
338 : RootFinderDataType &RootFinderData, // Data used by root finding algorithm
339 : Real64 const X, // X value of current iterate
340 : Real64 const Y, // Y value of current iterate
341 : bool &IsDoneFlag // If TRUE indicates that the iteration should be stopped
342 : )
343 : {
344 :
345 : // SUBROUTINE INFORMATION:
346 : // AUTHOR Dimitri Curtil (LBNL)
347 : // DATE WRITTEN March 2006
348 : // MODIFIED
349 : // RE-ENGINEERED na
350 :
351 : // PURPOSE OF THIS SUBROUTINE:
352 : // This subroutine is the workhorse of the root finder framework.
353 : // It should be invoked with every new candidate (X,Y) until
354 : // convergence is achieved or abnormal termination is detected.
355 : // The subroutine performs the following tasks:
356 : // - it checks for convergence
357 : // - it updates the internal data with the current iterate (X,Y)
358 : // - it computes a new root candidate if not converged yet.
359 : // Sets IsDoneFlag to FALSE when iteration should continue with the new candidate value.
360 : // Sets IsDoneFlag to TRUE when iteration should be stopped because convergence has been achieved
361 : // or because of a fatal error.
362 : // Note that only upon normal termination (iStatusOK<...> codes)
363 : // will the XCandidate value contain the root.
364 : // If the root has not been located yet the XCandidate value contains
365 : // the next candidate root to try.
366 : // Status IsDoneFlag XCandidate
367 : // ========================================================================
368 : // iStatusErrorRange TRUE na
369 : // iStatusErrorSingular TRUE na
370 : // iStatusErrorSlope TRUE na
371 : // iStatusErrorBracket TRUE na
372 : // ------------------------------------------------------------------------
373 : // iStatusOKMin TRUE MinPoint%X
374 : // iStatusOKMax TRUE MaxPoint%X
375 : // iStatusOK TRUE X
376 : // iStatusOKRoundOff TRUE X
377 : // ------------------------------------------------------------------------
378 : // iStatusNone FALSE AdvanceRootFinder()
379 : // iStatusWarningNonMonotonic FALSE AdvanceRootFinder()
380 : // iStatusWarningSingular FALSE AdvanceRootFinder()
381 :
382 : // METHODOLOGY EMPLOYED:
383 : // The methodology reflects the same approach implemented in the subroutine CalcSimpleController()
384 : // whereby the iteration was exited by checking the following conditions in this specified
385 : // sequence:
386 : // 1. Check for singular function so that YMin /= YMax
387 : // 2. Check for slope condition for the residual function
388 : // - increasing: YMin < YMax
389 : // - decreasing: YMin > YMax
390 : // 3. Check for min constraint
391 : // - increasing: YMin <= 0
392 : // - decreasing: YMin >= 0
393 : // 4. Check for max constraint
394 : // - increasing: YMax > 0
395 : // - decreasing: YMax < 0
396 : // 5. Check unconstrained convergence
397 : // Note that the slope condition was not explicitly checked in the original implementation
398 : // in CalcSimpleController().
399 : // Finally, we also check the X increments between successive iterations to detect possible
400 : // cases whereby the allowed precision in the X space limits the precision attainable in
401 : // the Y space. This check is implemented in:
402 : // - CheckIncrementRoundOff()
403 : // - CheckBracketRoundOff()
404 :
405 : // Reset status flag
406 19745373 : RootFinderData.StatusFlag = RootFinderStatus::None;
407 :
408 : // Check that MinPoint%X <= X <= MaxPoint%X
409 19745373 : if (!CheckMinMaxRange(RootFinderData, X)) {
410 0 : RootFinderData.StatusFlag = RootFinderStatus::ErrorRange;
411 :
412 : // Fatal error: No need to continue iterating
413 0 : IsDoneFlag = true;
414 0 : return;
415 : }
416 :
417 : // Update min/max support points with current iterate
418 19745373 : UpdateMinMax(RootFinderData, X, Y);
419 :
420 : //----------------------------------------------------------------------------
421 : // Check "global" singularity and bad slope conditions between min and
422 : // max points
423 : //----------------------------------------------------------------------------
424 :
425 : // NOTE: Performed before checking min and max constraints to mimic original implementation
426 : // in ManagerControllers()
427 19745373 : if (RootFinderData.MinPoint.DefinedFlag && RootFinderData.MaxPoint.DefinedFlag) {
428 :
429 : // Check that min and max points are distinct
430 2187848 : if (RootFinderData.MinPoint.X == RootFinderData.MaxPoint.X) {
431 440616 : RootFinderData.StatusFlag = RootFinderStatus::OKMin;
432 440616 : RootFinderData.XCandidate = RootFinderData.MinPoint.X;
433 :
434 : // Solution found: No need to continue iterating
435 440616 : IsDoneFlag = true;
436 440616 : return;
437 : }
438 :
439 1747232 : if (RootFinderData.MinPoint.DefinedFlag) {
440 1747232 : if (CheckMinConstraint(state, RootFinderData)) {
441 0 : RootFinderData.StatusFlag = RootFinderStatus::OKMin;
442 0 : RootFinderData.XCandidate = RootFinderData.MinPoint.X;
443 :
444 : // Solution found: No need to continue iterating
445 0 : IsDoneFlag = true;
446 0 : return;
447 : }
448 : }
449 :
450 : // Check singularity condition between min and max points
451 1747232 : if (!CheckNonSingularity(RootFinderData)) {
452 392731 : RootFinderData.StatusFlag = RootFinderStatus::ErrorSingular;
453 :
454 : // Fatal error: No need to continue iterating
455 392731 : IsDoneFlag = true;
456 392731 : return;
457 : }
458 :
459 : // Check slope condition between min and max points
460 1354501 : if (!CheckSlope(state, RootFinderData)) {
461 653 : RootFinderData.StatusFlag = RootFinderStatus::ErrorSlope;
462 :
463 : // Fatal error: No need to continue iterating
464 653 : IsDoneFlag = true;
465 653 : return;
466 : }
467 : }
468 :
469 : //----------------------------------------------------------------------------
470 : // Check that F(X) is not min or max constrained
471 : //----------------------------------------------------------------------------
472 :
473 : // Check min constraint before max constraint to mimic original implementation
474 : // in ManagerControllers()
475 18911373 : if (RootFinderData.MinPoint.DefinedFlag) {
476 18894659 : if (CheckMinConstraint(state, RootFinderData)) {
477 3621715 : RootFinderData.StatusFlag = RootFinderStatus::OKMin;
478 3621715 : RootFinderData.XCandidate = RootFinderData.MinPoint.X;
479 :
480 : // Solution found: No need to continue iterating
481 3621715 : IsDoneFlag = true;
482 3621715 : return;
483 : }
484 : }
485 :
486 : // Min point should always be evaluated first to ensure that we return with the min
487 : // consrained solution even in cases where the residual function has inconsistent slope.
488 : // TODO: Force to evaluate min point before exiting with max constrained solution
489 : // in order to be able to detect singularity and bad slope conditions.
490 15289658 : if (RootFinderData.MaxPoint.DefinedFlag) {
491 1356985 : if (CheckMaxConstraint(state, RootFinderData)) {
492 :
493 481208 : RootFinderData.StatusFlag = RootFinderStatus::OKMax;
494 481208 : RootFinderData.XCandidate = RootFinderData.MaxPoint.X;
495 :
496 : // Solution found: No need to continue iterating
497 481208 : IsDoneFlag = true;
498 481208 : return;
499 : }
500 : }
501 :
502 : //----------------------------------------------------------------------------
503 : // Check convergence of current iterate
504 : //----------------------------------------------------------------------------
505 :
506 : // Check unconstrained convergence after we are sure that the candidate X value lies
507 : // within the allowed min/max range
508 14808450 : if (CheckRootFinderConvergence(RootFinderData, Y)) {
509 3981728 : RootFinderData.StatusFlag = RootFinderStatus::OK;
510 3981728 : RootFinderData.XCandidate = X;
511 :
512 : // Update root finder internal data with current iterate (X,Y)
513 3981728 : UpdateRootFinder(state, RootFinderData, X, Y);
514 :
515 : // Solution found: No need to continue iterating
516 3981728 : IsDoneFlag = true;
517 3981728 : return;
518 : }
519 :
520 : // Check last as this was not done in the original implementation
521 : // Essentially we stop the iteration if:
522 : // - the distance between the lower and upper bounds is smaller than the user-specified
523 : // tolerance for the X variables. (USING brackets from previous iteration)
524 10826722 : if (CheckBracketRoundOff(RootFinderData)) {
525 3632 : RootFinderData.StatusFlag = RootFinderStatus::OKRoundOff;
526 :
527 : // Solution found: No need to continue iterating
528 3632 : IsDoneFlag = true;
529 3632 : return;
530 : }
531 :
532 : //----------------------------------------------------------------------------
533 : // If the current iterate lies within the current lower and upper brackets,
534 : // proceed with normal processing to identify the next root candidate:
535 : // - update lower/upper bracket with current iterate
536 : // - update history
537 : // - update increments across successive iterations
538 : // - update current point
539 : // - compute next candidate (see AdvanceRootFinder() ).
540 : //----------------------------------------------------------------------------
541 :
542 : // Check that current iterate is within the current lower and upper points
543 10823090 : if (!CheckLowerUpperBracket(RootFinderData, X)) {
544 0 : RootFinderData.StatusFlag = RootFinderStatus::ErrorBracket;
545 :
546 : // Fatal error: No need to continue iterating
547 0 : IsDoneFlag = true;
548 0 : return;
549 : }
550 :
551 : // Update root finder internal data with current iterate (X,Y)
552 10823090 : UpdateRootFinder(state, RootFinderData, X, Y);
553 :
554 : // Compute new root candidate and store value in in RootFinderData%XCandidate
555 : // - First attempt to bracket root within lower and upper points
556 : // - Then use whatever requested solution method in SetupRootFinder() to
557 : // compute the next candidate.
558 10823090 : AdvanceRootFinder(state, RootFinderData);
559 :
560 : // Indicates that we should continue iterating with new candidate
561 10823090 : IsDoneFlag = false;
562 : }
563 :
564 0 : RootFinderStatus CheckInternalConsistency(EnergyPlusData &state, RootFinderDataType const &RootFinderData) // Data used by root finding algorithm
565 : {
566 : // FUNCTION INFORMATION:
567 : // AUTHOR Dimitri Curtil (LBNL)
568 : // DATE WRITTEN March 2006
569 : // MODIFIED
570 : // RE-ENGINEERED na
571 :
572 : // PURPOSE OF THIS FUNCTION:
573 : // This function checks whether the lower and upper points (if defined)
574 : // determine a consistent interval bracketting the root.
575 : // Returns the status code accordingly.
576 : // This function does not modify the argument RooFinderData.
577 : // Only used internally for debugging.
578 :
579 : // Return value
580 : RootFinderStatus CheckInternalConsistency;
581 :
582 : // Default initialization
583 0 : CheckInternalConsistency = RootFinderStatus::None;
584 :
585 : // Internal consistency check involving both support points
586 0 : if (RootFinderData.LowerPoint.DefinedFlag && RootFinderData.UpperPoint.DefinedFlag) {
587 :
588 : // Check that the existing lower and upper points do bracket the root
589 0 : if (RootFinderData.LowerPoint.X > RootFinderData.UpperPoint.X) {
590 0 : CheckInternalConsistency = RootFinderStatus::ErrorRange;
591 0 : return CheckInternalConsistency;
592 : }
593 :
594 : // Check for non-monotonicity between the existing lower and upper points
595 0 : switch (RootFinderData.Controls.SlopeType) {
596 0 : case DataRootFinder::Slope::Increasing: {
597 : // Y-value of lower point must be strictly smaller than Y-value of upper point
598 0 : if (RootFinderData.LowerPoint.Y > RootFinderData.UpperPoint.Y) {
599 0 : CheckInternalConsistency = RootFinderStatus::WarningNonMonotonic;
600 0 : return CheckInternalConsistency;
601 : }
602 0 : } break;
603 0 : case DataRootFinder::Slope::Decreasing: {
604 : // Y-value of lower point must be strictly larger than Y-value of upper point
605 0 : if (RootFinderData.LowerPoint.Y < RootFinderData.UpperPoint.Y) {
606 0 : CheckInternalConsistency = RootFinderStatus::WarningNonMonotonic;
607 0 : return CheckInternalConsistency;
608 : }
609 0 : } break;
610 0 : default: {
611 : // Should never happen
612 0 : ShowSevereError(state, "CheckInternalConsistency: Invalid function slope specification. Valid choices are:");
613 0 : ShowContinueError(state, format("CheckInternalConsistency: Slope::Increasing={}", DataRootFinder::Slope::Increasing));
614 0 : ShowContinueError(state, format("CheckInternalConsistency: Slope::Decreasing={}", DataRootFinder::Slope::Decreasing));
615 0 : ShowFatalError(state, "CheckInternalConsistency: Preceding error causes program termination.");
616 0 : } break;
617 : }
618 :
619 : // Check for in singularity with respect to the existing lower and upper points
620 : // Only check if the lower and upper points are distinct!
621 0 : if (RootFinderData.UpperPoint.X > RootFinderData.LowerPoint.X) {
622 0 : if (RootFinderData.UpperPoint.Y == RootFinderData.LowerPoint.Y) {
623 0 : CheckInternalConsistency = RootFinderStatus::ErrorSingular;
624 0 : return CheckInternalConsistency;
625 : }
626 : }
627 : }
628 :
629 : // Check min constraint for min point if already defined
630 0 : if (RootFinderData.MinPoint.DefinedFlag) {
631 0 : switch (RootFinderData.Controls.SlopeType) {
632 0 : case DataRootFinder::Slope::Increasing: {
633 0 : if (RootFinderData.MinPoint.Y >= 0.0) {
634 0 : CheckInternalConsistency = RootFinderStatus::OKMin;
635 0 : return CheckInternalConsistency;
636 : }
637 0 : } break;
638 0 : case DataRootFinder::Slope::Decreasing: {
639 0 : if (RootFinderData.MinPoint.Y <= 0.0) {
640 0 : CheckInternalConsistency = RootFinderStatus::OKMin;
641 0 : return CheckInternalConsistency;
642 : }
643 0 : } break;
644 0 : default: {
645 : // Should never happen
646 0 : ShowSevereError(state, "CheckInternalConsistency: Invalid function slope specification. Valid choices are:");
647 0 : ShowContinueError(state, format("CheckInternalConsistency: Slope::Increasing={}", DataRootFinder::Slope::Increasing));
648 0 : ShowContinueError(state, format("CheckInternalConsistency: Slope::Decreasing={}", DataRootFinder::Slope::Decreasing));
649 0 : ShowFatalError(state, "CheckInternalConsistency: Preceding error causes program termination.");
650 0 : } break;
651 : }
652 : }
653 :
654 : // Check max constraint for max point if already defined
655 0 : if (RootFinderData.MaxPoint.DefinedFlag) {
656 0 : switch (RootFinderData.Controls.SlopeType) {
657 0 : case DataRootFinder::Slope::Increasing: {
658 0 : if (RootFinderData.MaxPoint.Y <= 0.0) {
659 0 : CheckInternalConsistency = RootFinderStatus::OKMax;
660 0 : return CheckInternalConsistency;
661 : }
662 0 : } break;
663 0 : case DataRootFinder::Slope::Decreasing: {
664 0 : if (RootFinderData.MaxPoint.Y >= 0.0) {
665 0 : CheckInternalConsistency = RootFinderStatus::OKMax;
666 0 : return CheckInternalConsistency;
667 : }
668 0 : } break;
669 0 : default: {
670 : // Should never happen
671 0 : ShowSevereError(state, "CheckInternalConsistency: Invalid function slope specification. Valid choices are:");
672 0 : ShowContinueError(state, format("CheckInternalConsistency: Slope::Increasing={}", DataRootFinder::Slope::Increasing));
673 0 : ShowContinueError(state, format("CheckInternalConsistency: Slope::Decreasing={}", DataRootFinder::Slope::Decreasing));
674 0 : ShowFatalError(state, "CheckInternalConsistency: Preceding error causes program termination.");
675 0 : } break;
676 : }
677 : }
678 :
679 0 : return CheckInternalConsistency;
680 : }
681 :
682 17589197 : bool CheckRootFinderCandidate(RootFinderDataType const &RootFinderData, // Data used by root finding algorithm
683 : Real64 const X // X value for current iterate
684 : )
685 : {
686 : // FUNCTION INFORMATION:
687 : // AUTHOR Dimitri Curtil (LBNL)
688 : // DATE WRITTEN February 2006
689 : // MODIFIED
690 : // RE-ENGINEERED na
691 :
692 : // PURPOSE OF THIS FUNCTION:
693 : // This function checks whether the root candidate X lies within the specified
694 : // min/max limits as well as within the current lower and upper brackets (if defined).
695 : // Returns TRUE if X value is a valid root candidate.
696 : // Returns FALSE otherwise.
697 :
698 : // Return value
699 : bool CheckRootFinderCandidate;
700 :
701 17589197 : if (CheckMinMaxRange(RootFinderData, X) && CheckLowerUpperBracket(RootFinderData, X)) {
702 17585761 : CheckRootFinderCandidate = true;
703 : } else {
704 3436 : CheckRootFinderCandidate = false;
705 : }
706 :
707 17589197 : return CheckRootFinderCandidate;
708 : }
709 :
710 37334570 : bool CheckMinMaxRange(RootFinderDataType const &RootFinderData, // Data used by root finding algorithm
711 : Real64 const X // X value for current iterate
712 : )
713 : {
714 : // FUNCTION INFORMATION:
715 : // AUTHOR Dimitri Curtil (LBNL)
716 : // DATE WRITTEN February 2006
717 : // MODIFIED Brent Griffith (NREL) added DefinedFlag traps
718 : // RE-ENGINEERED na
719 :
720 : // PURPOSE OF THIS FUNCTION:
721 : // This function checks whether the current iterate X lies within the specified min/max limits
722 : // or not.
723 : // Returns TRUE if current iterate satisfies min/max constraints.
724 : // Returns FALSE if current iterate is out-of-range.
725 :
726 : // Return value
727 : bool CheckMinMaxRange;
728 :
729 37334570 : if (RootFinderData.MinPoint.DefinedFlag) {
730 19597368 : if (X < RootFinderData.MinPoint.X) {
731 3227 : CheckMinMaxRange = false;
732 3227 : return CheckMinMaxRange;
733 : }
734 : }
735 :
736 37331343 : if (RootFinderData.MaxPoint.DefinedFlag) {
737 1304487 : if (X > RootFinderData.MaxPoint.X) {
738 2 : CheckMinMaxRange = false;
739 2 : return CheckMinMaxRange;
740 : }
741 : }
742 :
743 37331341 : CheckMinMaxRange = true;
744 :
745 37331341 : return CheckMinMaxRange;
746 : }
747 :
748 28409058 : bool CheckLowerUpperBracket(RootFinderDataType const &RootFinderData, // Data used by root finding algorithm
749 : Real64 const X // X value for current iterate
750 : )
751 : {
752 : // FUNCTION INFORMATION:
753 : // AUTHOR Dimitri Curtil (LBNL)
754 : // DATE WRITTEN February 2006
755 : // MODIFIED Brent Griffith, March 2010, changed to LowerPoint%X <= X <= UpperPoint%X
756 : // RE-ENGINEERED na
757 :
758 : // PURPOSE OF THIS FUNCTION:
759 : // This function checks whether the current iterate X lies within the current lower
760 : // and upper points or not (if defined):
761 : // LowerPoint%X < X < UpperPoint%X
762 : // Returns TRUE if current iterate lies within the lower/upper bracket.
763 : // Returns FALSE otherwise.
764 :
765 : // Return value
766 : bool CheckLowerUpperBracket;
767 :
768 28409058 : if (RootFinderData.LowerPoint.DefinedFlag) {
769 15543595 : if (X < RootFinderData.LowerPoint.X) {
770 89 : CheckLowerUpperBracket = false;
771 89 : return CheckLowerUpperBracket;
772 : }
773 : }
774 :
775 28408969 : if (RootFinderData.UpperPoint.DefinedFlag) {
776 7644414 : if (X > RootFinderData.UpperPoint.X) {
777 118 : CheckLowerUpperBracket = false;
778 118 : return CheckLowerUpperBracket;
779 : }
780 : }
781 :
782 28408851 : CheckLowerUpperBracket = true;
783 :
784 28408851 : return CheckLowerUpperBracket;
785 : }
786 :
787 1354501 : bool CheckSlope(EnergyPlusData &state, RootFinderDataType const &RootFinderData) // Data used by root finding algorithm
788 : {
789 : // FUNCTION INFORMATION:
790 : // AUTHOR Dimitri Curtil (LBNL)
791 : // DATE WRITTEN February 2006
792 : // MODIFIED
793 : // RE-ENGINEERED na
794 :
795 : // PURPOSE OF THIS FUNCTION:
796 : // This function checks whether the current iterate (X,Y) satisfies the slope
797 : // requirement between the min and max points.
798 : // Returns FALSE if the slope requirement is NOT satisfied.
799 : // Returns TRUE if the slope requirement is satisfied.
800 : // PRECONDITION:
801 : // - Function assumes that both the min and max points are defined.
802 : // POSTCONDITION:
803 : // - RootFinderData is NOT changed by this function.
804 :
805 : // Return value
806 : bool CheckSlope;
807 :
808 : // Check that the slope requirement is respected at the min and max points
809 : // Note that the singularity check takes care of RootFinderData%MinPoint%Y == RootFinderData%MaxPoint%Y
810 : // therefore we use strict comparison operators < and >.
811 1354501 : switch (RootFinderData.Controls.SlopeType) {
812 461679 : case DataRootFinder::Slope::Increasing: {
813 461679 : if (RootFinderData.MinPoint.Y < RootFinderData.MaxPoint.Y) {
814 461639 : CheckSlope = true;
815 461639 : return CheckSlope;
816 : }
817 40 : } break;
818 892822 : case DataRootFinder::Slope::Decreasing: {
819 892822 : if (RootFinderData.MinPoint.Y > RootFinderData.MaxPoint.Y) {
820 892209 : CheckSlope = true;
821 892209 : return CheckSlope;
822 : }
823 613 : } break;
824 0 : default: {
825 : // Should never happen
826 0 : ShowSevereError(state, "CheckSlope: Invalid function slope specification. Valid choices are:");
827 0 : ShowContinueError(state, format("CheckSlope: Slope::Increasing={}", DataRootFinder::Slope::Increasing));
828 0 : ShowContinueError(state, format("CheckSlope: Slope::Decreasing={}", DataRootFinder::Slope::Decreasing));
829 0 : ShowFatalError(state, "CheckSlope: Preceding error causes program termination.");
830 0 : } break;
831 : }
832 :
833 653 : CheckSlope = false;
834 :
835 653 : return CheckSlope;
836 : }
837 :
838 1747232 : bool CheckNonSingularity(RootFinderDataType const &RootFinderData) // Data used by root finding algorithm
839 : {
840 : // FUNCTION INFORMATION:
841 : // AUTHOR Dimitri Curtil (LBNL)
842 : // DATE WRITTEN February 2006
843 : // MODIFIED
844 : // RE-ENGINEERED na
845 :
846 : // PURPOSE OF THIS FUNCTION:
847 : // This function checks whether the min and max points define a locally, singular
848 : // equation system. In a 1-dimensional system, "singularity" is detected if the
849 : // min and max points have the same y-values, thereby producing a zero slope
850 : // across the min/max range.
851 : // Returns TRUE if the function satisfies the non-singularity condition.
852 : // Returns FALSE otherwise (i.e., F(X) essentially displays a zero slope
853 : // between the min and max points) .
854 : // PRECONDITION:
855 : // - Function assumes that both the min and max points are defined.
856 : // POSTCONDITION:
857 : // - RootFinderData is NOT changed by this function.
858 :
859 : // Return value
860 : bool CheckNonSingularity;
861 :
862 : // FUNCTION PARAMETER DEFINITIONS:
863 : // Safety factor used to detect a singular residual function between the min and max
864 : // points.
865 : // NOTE: Requesting exactly the same value is obtained by setting SafetyFactor = 0.0
866 1747232 : Real64 constexpr SafetyFactor(0.1);
867 :
868 : // FUNCTION LOCAL VARIABLE DECLARATIONS:
869 : Real64 DeltaY; // Difference between min and max Y-values
870 : Real64 ATolY; // Absolute tolerance used to detected equal min and max Y-values
871 :
872 : // Added this check based on an absolute tolerance test for y values to avoid incorrectly detecting
873 : // functions with bad slope due to numerical noise.
874 : // Typically, this takes care of situations where the controlled equipment in ManageControllers()
875 : // would be misdiagnosed as displaying the "wrong slope" instead of being treated as "singular"
876 : // (i.e. in inactive mode).
877 1747232 : DeltaY = std::abs(RootFinderData.MinPoint.Y - RootFinderData.MaxPoint.Y);
878 1747232 : ATolY = SafetyFactor * RootFinderData.Controls.ATolY;
879 :
880 1747232 : if (std::abs(DeltaY) <= ATolY) {
881 392731 : CheckNonSingularity = false;
882 : } else {
883 1354501 : CheckNonSingularity = true;
884 : }
885 :
886 1747232 : return CheckNonSingularity;
887 : }
888 :
889 20641891 : bool CheckMinConstraint(EnergyPlusData &state, RootFinderDataType const &RootFinderData) // Data used by root finding algorithm
890 : {
891 : // FUNCTION INFORMATION:
892 : // AUTHOR Dimitri Curtil (LBNL)
893 : // DATE WRITTEN February 2006
894 : // MODIFIED
895 : // RE-ENGINEERED na
896 :
897 : // PURPOSE OF THIS FUNCTION:
898 : // This function checks whether the min point satisfies the min constraint
899 : // condition or not.
900 : // PRECONDITION:
901 : // - Function assumes that the min point is defined.
902 : // POSTCONDITION:
903 : // - RootFinderData is NOT changed by this function.
904 :
905 : // Return value
906 : bool CheckMinConstraint;
907 :
908 20641891 : switch (RootFinderData.Controls.SlopeType) {
909 6693551 : case DataRootFinder::Slope::Increasing: {
910 6693551 : if (RootFinderData.MinPoint.Y >= 0.0) {
911 2846801 : CheckMinConstraint = true;
912 2846801 : return CheckMinConstraint;
913 : }
914 3846750 : } break;
915 13948340 : case DataRootFinder::Slope::Decreasing: {
916 13948340 : if (RootFinderData.MinPoint.Y <= 0.0) {
917 774914 : CheckMinConstraint = true;
918 774914 : return CheckMinConstraint;
919 : }
920 13173426 : } break;
921 0 : default: {
922 : // Should never happen
923 0 : ShowSevereError(state, "CheckMinConstraint: Invalid function slope specification. Valid choices are:");
924 0 : ShowContinueError(state, format("CheckMinConstraint: Slope::Increasing={}", DataRootFinder::Slope::Increasing));
925 0 : ShowContinueError(state, format("CheckMinConstraint: Slope::Decreasing={}", DataRootFinder::Slope::Decreasing));
926 0 : ShowFatalError(state, "CheckMinConstraint: Preceding error causes program termination.");
927 0 : } break;
928 : }
929 :
930 17020176 : CheckMinConstraint = false;
931 :
932 17020176 : return CheckMinConstraint;
933 : }
934 :
935 1356985 : bool CheckMaxConstraint(EnergyPlusData &state, RootFinderDataType const &RootFinderData) // Data used by root finding algorithm
936 : {
937 : // FUNCTION INFORMATION:
938 : // AUTHOR Dimitri Curtil (LBNL)
939 : // DATE WRITTEN February 2006
940 : // MODIFIED
941 : // RE-ENGINEERED na
942 :
943 : // PURPOSE OF THIS FUNCTION:
944 : // This function checks whether the max point satisfies the max constraint
945 : // condition or not.
946 : // PRECONDITION:
947 : // - Function assumes that the max point is defined.
948 : // POSTCONDITION:
949 : // - RootFinderData is NOT changed by this function.
950 :
951 : // Return value
952 : bool CheckMaxConstraint;
953 :
954 : // Check for max constrained convergence with respect to the new iterate (X,Y)
955 1356985 : switch (RootFinderData.Controls.SlopeType) {
956 464776 : case DataRootFinder::Slope::Increasing: {
957 464776 : if (RootFinderData.MaxPoint.Y <= 0.0) {
958 147261 : CheckMaxConstraint = true;
959 147261 : return CheckMaxConstraint;
960 : }
961 317515 : } break;
962 892209 : case DataRootFinder::Slope::Decreasing: {
963 892209 : if (RootFinderData.MaxPoint.Y >= 0.0) {
964 333947 : CheckMaxConstraint = true;
965 333947 : return CheckMaxConstraint;
966 : }
967 558262 : } break;
968 0 : default: {
969 : // Should never happen
970 0 : ShowSevereError(state, "CheckMaxConstraint: Invalid function slope specification. Valid choices are:");
971 0 : ShowContinueError(state, format("CheckMaxConstraint: Slope::Increasing={}", DataRootFinder::Slope::Increasing));
972 0 : ShowContinueError(state, format("CheckMaxConstraint: Slope::Decreasing={}", DataRootFinder::Slope::Decreasing));
973 0 : ShowFatalError(state, "CheckMaxConstraint: Preceding error causes program termination.");
974 0 : } break;
975 : }
976 :
977 875777 : CheckMaxConstraint = false;
978 :
979 875777 : return CheckMaxConstraint;
980 : }
981 :
982 21662697 : bool CheckRootFinderConvergence(RootFinderDataType const &RootFinderData, // Data used by root finding algorithm
983 : Real64 const Y // Y value for current iterate
984 : )
985 : {
986 : // FUNCTION INFORMATION:
987 : // AUTHOR Dimitri Curtil (LBNL)
988 : // DATE WRITTEN February 2006
989 : // MODIFIED
990 : // RE-ENGINEERED na
991 :
992 : // PURPOSE OF THIS FUNCTION:
993 : // This function checks whether the current iterate (X,Y) satisfies the
994 : // unconstrained convergence criterion or not.
995 :
996 : // Return value
997 : bool CheckRootFinderConvergence;
998 :
999 : // Check for unconstrained convergence
1000 21662697 : if (std::abs(Y) <= RootFinderData.Controls.ATolY) {
1001 8599979 : CheckRootFinderConvergence = true;
1002 8599979 : return CheckRootFinderConvergence;
1003 : }
1004 :
1005 13062718 : CheckRootFinderConvergence = false;
1006 :
1007 13062718 : return CheckRootFinderConvergence;
1008 : }
1009 :
1010 10826722 : bool CheckBracketRoundOff(RootFinderDataType const &RootFinderData) // Data used by root finding algorithm
1011 : {
1012 : // FUNCTION INFORMATION:
1013 : // AUTHOR Dimitri Curtil (LBNL)
1014 : // DATE WRITTEN February 2006
1015 : // MODIFIED
1016 : // RE-ENGINEERED na
1017 :
1018 : // PURPOSE OF THIS FUNCTION:
1019 : // This function checks whether the current lower and upper brackets satisfies
1020 : // the round-off criterion or not.
1021 :
1022 : // Return value
1023 : bool CheckBracketRoundOff;
1024 :
1025 : // FUNCTION LOCAL VARIABLE DECLARATIONS:
1026 : Real64 DeltaUL; // Distance between lower and upper points
1027 : Real64 TypUL; // Typical value for values lying within lower/upper interval
1028 : Real64 TolUL; // Tolerance to satisfy for lower-upper distance
1029 :
1030 : // Check for round-off error in Lower/Upper interval
1031 10826722 : if (RootFinderData.LowerPoint.DefinedFlag && RootFinderData.UpperPoint.DefinedFlag) {
1032 3644361 : DeltaUL = RootFinderData.UpperPoint.X - RootFinderData.LowerPoint.X;
1033 3644361 : TypUL = (std::abs(RootFinderData.UpperPoint.X) + std::abs(RootFinderData.LowerPoint.X)) / 2.0;
1034 3644361 : TolUL = RootFinderData.Controls.TolX * std::abs(TypUL) + RootFinderData.Controls.ATolX;
1035 :
1036 : // Halve tolerance to reflect the fact that solution can be anywhere between the lower and upper points.
1037 3644361 : if (std::abs(DeltaUL) <= 0.5 * std::abs(TolUL)) {
1038 3632 : CheckBracketRoundOff = true;
1039 3632 : return CheckBracketRoundOff;
1040 : }
1041 : }
1042 :
1043 10823090 : CheckBracketRoundOff = false;
1044 :
1045 10823090 : return CheckBracketRoundOff;
1046 : }
1047 :
1048 19745373 : void UpdateMinMax(RootFinderDataType &RootFinderData, // Data used by root finding algorithm
1049 : Real64 const X, // X value for current iterate
1050 : Real64 const Y // Y value for current iterate, F(X)=Y
1051 : )
1052 : {
1053 : // SUBROUTINE INFORMATION:
1054 : // AUTHOR Dimitri Curtil (LBNL)
1055 : // DATE WRITTEN February 2006
1056 : // MODIFIED
1057 : // RE-ENGINEERED na
1058 :
1059 : // PURPOSE OF THIS SUBROUTINE:
1060 : // This subroutine updates the min/max support points in the root finder data.
1061 :
1062 : // METHODOLOGY EMPLOYED:
1063 : // PRECONDITION:
1064 : // na
1065 : // POSTCONDITION:
1066 : // - RootFinderData%MinPoint possibly updated
1067 : // - RootFinderData%MaxPoint possibly updated
1068 :
1069 : // Update min support point
1070 19745373 : if (X == RootFinderData.MinPoint.X) {
1071 8915016 : RootFinderData.MinPoint.Y = Y;
1072 8915016 : RootFinderData.MinPoint.DefinedFlag = true;
1073 : }
1074 :
1075 : // Update max support point
1076 19745373 : if (X == RootFinderData.MaxPoint.X) {
1077 1464189 : RootFinderData.MaxPoint.Y = Y;
1078 1464189 : RootFinderData.MaxPoint.DefinedFlag = true;
1079 : }
1080 19745373 : }
1081 :
1082 14804818 : void UpdateBracket(EnergyPlusData &state,
1083 : RootFinderDataType &RootFinderData, // Data used by root finding algorithm
1084 : Real64 const X, // X value for current iterate
1085 : Real64 const Y // Y value for current iterate, F(X)=Y
1086 : )
1087 : {
1088 : // SUBROUTINE INFORMATION:
1089 : // AUTHOR Dimitri Curtil (LBNL)
1090 : // DATE WRITTEN March 2006
1091 : // MODIFIED
1092 : // RE-ENGINEERED na
1093 :
1094 : // PURPOSE OF THIS SUBROUTINE:
1095 : // This subroutine updates the lower/upper support points in the root finder data
1096 : // with the current iterate (X,Y).
1097 :
1098 : // METHODOLOGY EMPLOYED:
1099 : // PRECONDITION:
1100 : // - The current iterate (X,Y) must satisfy:
1101 : // MinPoint%X <= LowerPoint%X < X < UpperPoint%X <= MaxPoint%X
1102 : // - RootFinderData%StatusFlag == iStatusNone
1103 : // POSTCONDITION:
1104 : // - RootFinderData%LowerPoint possibly updated
1105 : // - RootFinderData%UpperPoint possibly updated
1106 : // - RootFinderData%StatusFlag possibly updated with:
1107 : // - iStatusWarningNonMonotonic
1108 : // - iStatusWarningSingular
1109 :
1110 14804818 : switch (RootFinderData.Controls.SlopeType) {
1111 3186567 : case DataRootFinder::Slope::Increasing: {
1112 : // Update lower point
1113 3186567 : if (Y <= 0.0) {
1114 2235410 : if (!RootFinderData.LowerPoint.DefinedFlag) {
1115 1379608 : RootFinderData.LowerPoint.DefinedFlag = true;
1116 1379608 : RootFinderData.LowerPoint.X = X;
1117 1379608 : RootFinderData.LowerPoint.Y = Y;
1118 : } else {
1119 855802 : if (X >= RootFinderData.LowerPoint.X) {
1120 855802 : if (Y == RootFinderData.LowerPoint.Y) {
1121 1092 : RootFinderData.StatusFlag = RootFinderStatus::WarningSingular;
1122 854710 : } else if (Y < RootFinderData.LowerPoint.Y) {
1123 1081 : RootFinderData.StatusFlag = RootFinderStatus::WarningNonMonotonic;
1124 : }
1125 : // Update lower point with current iterate
1126 855802 : RootFinderData.LowerPoint.X = X;
1127 855802 : RootFinderData.LowerPoint.Y = Y;
1128 : } else {
1129 : // Should never happen if CheckLowerUpperBracket() is called before
1130 0 : ShowSevereError(state, "UpdateBracket: Current iterate is smaller than the lower bracket.");
1131 0 : ShowContinueError(state, format("UpdateBracket: X={:.15T}, Y={:.15T}", X, Y));
1132 0 : ShowContinueError(
1133 0 : state, format("UpdateBracket: XLower={:.15T}, YLower={:.15T}", RootFinderData.LowerPoint.X, RootFinderData.LowerPoint.Y));
1134 0 : ShowFatalError(state, "UpdateBracket: Preceding error causes program termination.");
1135 : }
1136 : }
1137 :
1138 : // Update upper point
1139 : } else {
1140 951157 : if (!RootFinderData.UpperPoint.DefinedFlag) {
1141 392561 : RootFinderData.UpperPoint.DefinedFlag = true;
1142 392561 : RootFinderData.UpperPoint.X = X;
1143 392561 : RootFinderData.UpperPoint.Y = Y;
1144 : } else {
1145 558596 : if (X <= RootFinderData.UpperPoint.X) {
1146 558596 : if (Y == RootFinderData.UpperPoint.Y) {
1147 0 : RootFinderData.StatusFlag = RootFinderStatus::WarningSingular;
1148 558596 : } else if (Y > RootFinderData.UpperPoint.Y) {
1149 721 : RootFinderData.StatusFlag = RootFinderStatus::WarningNonMonotonic;
1150 : }
1151 : // Update upper point with current iterate
1152 558596 : RootFinderData.UpperPoint.X = X;
1153 558596 : RootFinderData.UpperPoint.Y = Y;
1154 : } else {
1155 : // Should never happen if CheckLowerUpperBracket() is called before
1156 0 : ShowSevereError(state, "UpdateBracket: Current iterate is greater than the upper bracket.");
1157 0 : ShowContinueError(state, format("UpdateBracket: X={:.15T}, Y={:.15T}", X, Y));
1158 0 : ShowContinueError(
1159 0 : state, format("UpdateBracket: XUpper={:.15T}, YUpper={:.15T}", RootFinderData.UpperPoint.X, RootFinderData.UpperPoint.Y));
1160 0 : ShowFatalError(state, "UpdateBracket: Preceding error causes program termination.");
1161 : }
1162 : }
1163 : }
1164 : // Monotone, decreasing function
1165 3186567 : } break;
1166 11618251 : case DataRootFinder::Slope::Decreasing: {
1167 : // Update lower point
1168 11618251 : if (Y >= 0.0) {
1169 7902057 : if (!RootFinderData.LowerPoint.DefinedFlag) {
1170 3475860 : RootFinderData.LowerPoint.DefinedFlag = true;
1171 3475860 : RootFinderData.LowerPoint.X = X;
1172 3475860 : RootFinderData.LowerPoint.Y = Y;
1173 : } else {
1174 4426197 : if (X >= RootFinderData.LowerPoint.X) {
1175 4426197 : if (Y == RootFinderData.LowerPoint.Y) {
1176 13639 : RootFinderData.StatusFlag = RootFinderStatus::WarningSingular;
1177 4412558 : } else if (Y > RootFinderData.LowerPoint.Y) {
1178 917 : RootFinderData.StatusFlag = RootFinderStatus::WarningNonMonotonic;
1179 : }
1180 : // Update lower point with current iterate
1181 4426197 : RootFinderData.LowerPoint.X = X;
1182 4426197 : RootFinderData.LowerPoint.Y = Y;
1183 : } else {
1184 : // Should never happen if CheckLowerUpperBracket() is called before
1185 0 : ShowSevereError(state, "UpdateBracket: Current iterate is smaller than the lower bracket.");
1186 0 : ShowContinueError(state, format("UpdateBracket: X={:.15T}, Y={:.15T}", X, Y));
1187 0 : ShowContinueError(
1188 0 : state, format("UpdateBracket: XLower={:.15T}, YLower={:.15T}", RootFinderData.LowerPoint.X, RootFinderData.LowerPoint.Y));
1189 0 : ShowFatalError(state, "UpdateBracket: Preceding error causes program termination.");
1190 : }
1191 : }
1192 :
1193 : // Update upper point
1194 : } else {
1195 3716194 : if (!RootFinderData.UpperPoint.DefinedFlag) {
1196 1299973 : RootFinderData.UpperPoint.DefinedFlag = true;
1197 1299973 : RootFinderData.UpperPoint.X = X;
1198 1299973 : RootFinderData.UpperPoint.Y = Y;
1199 : } else {
1200 2416221 : if (X <= RootFinderData.UpperPoint.X) {
1201 2416221 : if (Y == RootFinderData.UpperPoint.Y) {
1202 114 : RootFinderData.StatusFlag = RootFinderStatus::WarningSingular;
1203 2416107 : } else if (Y < RootFinderData.UpperPoint.Y) {
1204 2601 : RootFinderData.StatusFlag = RootFinderStatus::WarningNonMonotonic;
1205 : }
1206 : // Update upper point with current iterate
1207 2416221 : RootFinderData.UpperPoint.X = X;
1208 2416221 : RootFinderData.UpperPoint.Y = Y;
1209 : } else {
1210 : // Should never happen if CheckLowerUpperBracket() is called before
1211 0 : ShowSevereError(state, "UpdateBracket: Current iterate is greater than the upper bracket.");
1212 0 : ShowContinueError(state, format("UpdateBracket: X={:.15T}, Y={:.15T}", X, Y));
1213 0 : ShowContinueError(
1214 0 : state, format("UpdateBracket: XUpper={:.15T}, YUpper={:.15T}", RootFinderData.UpperPoint.X, RootFinderData.UpperPoint.Y));
1215 0 : ShowFatalError(state, "UpdateBracket: Preceding error causes program termination.");
1216 : }
1217 : }
1218 : }
1219 11618251 : } break;
1220 0 : default: {
1221 : // Should never happen
1222 0 : ShowSevereError(state, "UpdateBracket: Invalid function slope specification. Valid choices are:");
1223 0 : ShowContinueError(state, format("UpdateBracket: Slope::Increasing={}", DataRootFinder::Slope::Increasing));
1224 0 : ShowContinueError(state, format("UpdateBracket: Slope::Decreasing={}", DataRootFinder::Slope::Decreasing));
1225 0 : ShowFatalError(state, "UpdateBracket: Preceding error causes program termination.");
1226 0 : } break;
1227 : }
1228 14804818 : }
1229 :
1230 14804818 : void UpdateHistory(RootFinderDataType &RootFinderData, // Data used by root finding algorithm
1231 : Real64 const X, // X value for current iterate
1232 : Real64 const Y // Y value for current iterate, F(X)=Y
1233 : )
1234 : {
1235 : // SUBROUTINE INFORMATION:
1236 : // AUTHOR Dimitri Curtil (LBNL)
1237 : // DATE WRITTEN February 2006
1238 : // MODIFIED
1239 : // RE-ENGINEERED na
1240 :
1241 : // PURPOSE OF THIS SUBROUTINE:
1242 : // This subroutine updates the min/max support points in the root finder data.
1243 :
1244 : // METHODOLOGY EMPLOYED:
1245 : // PRECONDITION:
1246 : // - The current iterate (X,Y) must be a valid iterate:
1247 : // MinPoint%X <= LowerPoint%X < X < UpperPoint%X <= MaxPoint%X
1248 : // POSTCONDITION:
1249 : // - RootFinderData%History(:) updated with last 3 best iterates
1250 :
1251 : // SUBROUTINE LOCAL VARIABLE DECLARATIONS:
1252 : int NumHistory;
1253 :
1254 : // Update history with best iterates so that:
1255 : // ABS(History(1)%Y) <= ABS(History(2)%Y) <= ABS(History(3)%Y)
1256 : // Note that the history points are sorted so that
1257 : // SIGN(History(1)%Y) = -SIGN(History(3)%Y)
1258 : // to ensure that the history points bracket the candidate root.
1259 59219272 : for (auto &e : RootFinderData.History) {
1260 44414454 : e.X = e.Y = 0.0;
1261 44414454 : e.DefinedFlag = false;
1262 : }
1263 :
1264 14804818 : NumHistory = 0;
1265 14804818 : if (RootFinderData.LowerPoint.DefinedFlag) {
1266 9943648 : ++NumHistory;
1267 9943648 : RootFinderData.History(NumHistory).DefinedFlag = RootFinderData.LowerPoint.DefinedFlag;
1268 9943648 : RootFinderData.History(NumHistory).X = RootFinderData.LowerPoint.X;
1269 9943648 : RootFinderData.History(NumHistory).Y = RootFinderData.LowerPoint.Y;
1270 : }
1271 14804818 : if (RootFinderData.UpperPoint.DefinedFlag) {
1272 5025425 : ++NumHistory;
1273 5025425 : RootFinderData.History(NumHistory).DefinedFlag = RootFinderData.UpperPoint.DefinedFlag;
1274 5025425 : RootFinderData.History(NumHistory).X = RootFinderData.UpperPoint.X;
1275 5025425 : RootFinderData.History(NumHistory).Y = RootFinderData.UpperPoint.Y;
1276 : }
1277 14804818 : ++NumHistory;
1278 14804818 : RootFinderData.History(NumHistory).DefinedFlag = true;
1279 14804818 : RootFinderData.History(NumHistory).X = X;
1280 14804818 : RootFinderData.History(NumHistory).Y = Y;
1281 :
1282 14804818 : RootFinderData.NumHistory = NumHistory;
1283 14804818 : SortHistory(NumHistory, RootFinderData.History);
1284 14804818 : }
1285 :
1286 14804818 : void UpdateRootFinder(EnergyPlusData &state,
1287 : RootFinderDataType &RootFinderData, // Data used by root finding algorithm
1288 : Real64 const X, // X value for current iterate
1289 : Real64 const Y // Y value for current iterate, F(X)=Y
1290 : )
1291 : {
1292 : // SUBROUTINE INFORMATION:
1293 : // AUTHOR Dimitri Curtil (LBNL)
1294 : // DATE WRITTEN February 2006
1295 : // MODIFIED
1296 : // RE-ENGINEERED na
1297 :
1298 : // PURPOSE OF THIS SUBROUTINE:
1299 : // This subroutine updates the root finder internal data to account for
1300 : // the current iterate (X,Y):
1301 : // - Lower / Upper support points
1302 : // - Increments for successive iterates
1303 : // - Convergence rate
1304 :
1305 : // METHODOLOGY EMPLOYED:
1306 : // PRECONDITION:
1307 : // - The current iterate (X,Y) must be a valid iterate:
1308 : // MinPoint%X <= LowerPoint%X < X < UpperPoint%X <= MaxPoint%X
1309 : // - Invoke UpdateRootFinder() only if:
1310 : // - CheckRootFinderCandidate() returned TRUE
1311 : // - CheckNonSingularity() returned TRUE
1312 : // - CheckSlope() returned TRUE
1313 : // POSTCONDITION:
1314 : // - RootFinderData%LowerPoint possibly updated
1315 : // - RootFinderData%UpperPoint possibly updated
1316 : // - RootFinderData%CurrentPoint updated with current iterate (X,Y)
1317 : // - RootFinderData%History(:) updated with last 3 best iterates
1318 : // - RootFinderData%Increment updated
1319 : // - RootFinderData%ConvergenceRate updated
1320 :
1321 : // Update history with best iterates so that:
1322 : // ABS(History(1)%Y) <= ABS(History(2)%Y) <= ABS(History(3)%Y)
1323 : // Note that we must update the history before updating the lower/upper points
1324 14804818 : UpdateHistory(RootFinderData, X, Y);
1325 :
1326 : // Update lower and upper points
1327 14804818 : UpdateBracket(state, RootFinderData, X, Y);
1328 :
1329 : // Update increments and convergence rate
1330 14804818 : if (RootFinderData.CurrentPoint.DefinedFlag) {
1331 9943648 : RootFinderData.Increment.DefinedFlag = true;
1332 9943648 : RootFinderData.Increment.X = X - RootFinderData.CurrentPoint.X;
1333 9943648 : RootFinderData.Increment.Y = Y - RootFinderData.CurrentPoint.Y;
1334 :
1335 9943648 : if (std::abs(RootFinderData.CurrentPoint.Y) > 0.0) {
1336 : // NOTE: Should be smaller than one for convergent process
1337 9943648 : RootFinderData.ConvergenceRate = std::abs(Y) / std::abs(RootFinderData.CurrentPoint.Y);
1338 : } else {
1339 : // NOTE: Should never happen
1340 0 : RootFinderData.ConvergenceRate = -1.0;
1341 : }
1342 : }
1343 :
1344 : // Finally update CurrentPoint (must be done last)
1345 14804818 : RootFinderData.CurrentPoint.DefinedFlag = true;
1346 14804818 : RootFinderData.CurrentPoint.X = X;
1347 14804818 : RootFinderData.CurrentPoint.Y = Y;
1348 14804818 : }
1349 :
1350 14804818 : void SortHistory(int const N, // Number of points to sort in history array
1351 : Array1D<PointType> &History // Array of PointType variables. At least N of them
1352 : )
1353 : {
1354 : // SUBROUTINE INFORMATION:
1355 : // AUTHOR Dimitri Curtil (LBNL)
1356 : // DATE WRITTEN March 2006
1357 : // MODIFIED
1358 : // RE-ENGINEERED na
1359 :
1360 : // PURPOSE OF THIS SUBROUTINE:
1361 : // This subroutine orders the N points in the history array in increasing
1362 : // order of ABS(Y) values.
1363 :
1364 : // SUBROUTINE LOCAL VARIABLE DECLARATIONS:
1365 : int I;
1366 : int J;
1367 : Real64 XTemp;
1368 : Real64 YTemp;
1369 :
1370 : // Nothing to do if only one point stored in history
1371 14804818 : if (N <= 1) {
1372 4861170 : return;
1373 : }
1374 :
1375 24912721 : for (I = 1; I <= N - 1; ++I) {
1376 34963571 : for (J = I + 1; J <= N; ++J) {
1377 19994498 : if (History(J).DefinedFlag) {
1378 : // Swap I and J elements
1379 19994498 : if (std::abs(History(J).Y) < std::abs(History(I).Y)) {
1380 18459350 : XTemp = History(I).X;
1381 18459350 : YTemp = History(I).Y;
1382 18459350 : History(I).X = History(J).X;
1383 18459350 : History(I).Y = History(J).Y;
1384 18459350 : History(J).X = XTemp;
1385 18459350 : History(J).Y = YTemp;
1386 : }
1387 : }
1388 : }
1389 : }
1390 : }
1391 :
1392 10823090 : void AdvanceRootFinder(EnergyPlusData &state, RootFinderDataType &RootFinderData) // Data used by root finding algorithm
1393 : {
1394 : // SUBROUTINE INFORMATION:
1395 : // AUTHOR Dimitri Curtil (LBNL)
1396 : // DATE WRITTEN February 2006
1397 : // MODIFIED
1398 : // RE-ENGINEERED na
1399 :
1400 : // PURPOSE OF THIS SUBROUTINE:
1401 : // This subroutine computes the next candidate value based on the information available so far.
1402 : // Stores new value into RootFinderData%XCandidate
1403 : // PRECONDITION:
1404 : // na
1405 : // POSTCONDITION:
1406 : // - LowerPoint%X < XCandidate < UpperPoint%X
1407 : // - RootFinderData%CurrentMethodType update with current solution method.
1408 :
1409 : // METHODOLOGY EMPLOYED:
1410 : // The subroutine first attempts to bracket the root within a lower and upper point.
1411 : // Once it is bracketed, then we use the specified solution methods (Bisection,
1412 : // False position, Secant and Brent) to compute the next candidate.
1413 :
1414 : // SUBROUTINE LOCAL VARIABLE DECLARATIONS:
1415 10823090 : auto &XNext = state.dataGeneral->XNext;
1416 :
1417 : //----------------------------------------------------------------------------
1418 : // First attempt to bracket root between a lower point and an upper point.
1419 : //----------------------------------------------------------------------------
1420 :
1421 : // Detect the lower bracket
1422 10823090 : if (!RootFinderData.LowerPoint.DefinedFlag) {
1423 1169 : RootFinderData.CurrentMethodType = RootFinderMethod::Bracket;
1424 : // If we have 2 points already, try to detect lower point using the Secant formula
1425 1169 : if (BracketRoot(RootFinderData, XNext)) {
1426 0 : RootFinderData.XCandidate = XNext;
1427 : } else {
1428 1169 : if (!RootFinderData.MinPoint.DefinedFlag) {
1429 1169 : RootFinderData.XCandidate = RootFinderData.MinPoint.X;
1430 : } else {
1431 : // Should never happen
1432 0 : ShowFatalError(state, "AdvanceRootFinder: Cannot find lower bracket.");
1433 : }
1434 : }
1435 :
1436 : // Detect the upper bracket
1437 10821921 : } else if (!RootFinderData.UpperPoint.DefinedFlag) {
1438 5792849 : RootFinderData.CurrentMethodType = RootFinderMethod::Bracket;
1439 : // If we have 2 points already, try to detect upper point using the Secant formula
1440 5792849 : if (BracketRoot(RootFinderData, XNext)) {
1441 1727905 : RootFinderData.XCandidate = XNext;
1442 : } else {
1443 4064944 : if (!RootFinderData.MaxPoint.DefinedFlag) {
1444 4064944 : RootFinderData.XCandidate = RootFinderData.MaxPoint.X;
1445 : } else {
1446 : // Should never happen
1447 0 : ShowFatalError(state, "AdvanceRootFinder: Cannot find upper bracket.");
1448 : }
1449 : }
1450 :
1451 : //----------------------------------------------------------------------------
1452 : // Root finding can start ...
1453 : // Assumptions:
1454 : // - the lower and upper support points are defined.
1455 : // - the increments are defined (at least 2 history points are available)
1456 : //----------------------------------------------------------------------------
1457 : } else {
1458 5029072 : switch (RootFinderData.StatusFlag) {
1459 0 : case RootFinderStatus::OKRoundOff: {
1460 : // Should never happen if we exit the root finder upon detecting round-off condition
1461 0 : RootFinderData.XCandidate = BisectionMethod(RootFinderData);
1462 0 : } break;
1463 14343 : case RootFinderStatus::WarningSingular:
1464 : case RootFinderStatus::WarningNonMonotonic: {
1465 : // Following local singularity or non-monotonicity warnings we attempt
1466 : // to recover with the false position method to avoid running into trouble
1467 : // because the latest iterate did nt produce any improvement compared to
1468 : // the previous lower and upper brackets.
1469 14343 : RootFinderData.XCandidate = FalsePositionMethod(RootFinderData);
1470 14343 : } break;
1471 5014729 : default: {
1472 : // Assuming that the root is bracketed between the lower and upper points,
1473 : // we execute the requested solution method to produce the next candidate value
1474 : // for the root.
1475 5014729 : switch (RootFinderData.Controls.MethodType) {
1476 0 : case RootFinderMethod::Bisection: {
1477 : // Bisection method (aka interval halving)
1478 0 : RootFinderData.XCandidate = BisectionMethod(RootFinderData);
1479 0 : } break;
1480 38642 : case RootFinderMethod::FalsePosition: {
1481 : // False position method (aka regula falsi)
1482 38642 : RootFinderData.XCandidate = FalsePositionMethod(RootFinderData);
1483 38642 : } break;
1484 0 : case RootFinderMethod::Secant: {
1485 : // Secant method
1486 0 : RootFinderData.XCandidate = SecantMethod(RootFinderData);
1487 0 : } break;
1488 4976087 : case RootFinderMethod::Brent: {
1489 : // Brent method
1490 4976087 : RootFinderData.XCandidate = BrentMethod(RootFinderData);
1491 4976087 : } break;
1492 0 : default: {
1493 0 : ShowSevereError(state, "AdvanceRootFinder: Invalid solution method specification. Valid choices are:");
1494 0 : ShowContinueError(state, format("AdvanceRootFinder: iMethodBisection={}", RootFinderMethod::Bisection));
1495 0 : ShowContinueError(state, format("AdvanceRootFinder: iMethodFalsePosition={}", RootFinderMethod::FalsePosition));
1496 0 : ShowContinueError(state, format("AdvanceRootFinder: iMethodSecant={}", RootFinderMethod::Secant));
1497 0 : ShowContinueError(state, format("AdvanceRootFinder: iMethodBrent={}", RootFinderMethod::Brent));
1498 0 : ShowFatalError(state, "AdvanceRootFinder: Preceding error causes program termination.");
1499 0 : } break;
1500 : }
1501 5014729 : } break;
1502 : }
1503 : }
1504 10823090 : }
1505 :
1506 5794018 : bool BracketRoot(RootFinderDataType const &RootFinderData, // Data used by root finding algorithm
1507 : Real64 &XNext // Next value
1508 : )
1509 : {
1510 : // FUNCTION INFORMATION:
1511 : // AUTHOR Dimitri Curtil (LBNL)
1512 : // DATE WRITTEN February 2006
1513 : // MODIFIED
1514 : // RE-ENGINEERED na
1515 :
1516 : // PURPOSE OF THIS FUNCTION:
1517 : // This function attempts to compute a new point that will bracket the root
1518 : // using the secant formula to take advantage of the slope between the last 2
1519 : // iterates.
1520 : // Returns TRUE if successfully computed a new bracket in XNext.
1521 : // Else returns FASLE and does not update the XNext argument.
1522 : // Should only be used while in braketing mode (iMethodBracket).
1523 : // When the lower and upper brackets are detected then the FUNCTION SecantMethod
1524 : // should be used instead.
1525 : // PRECONDITION:
1526 : // na
1527 : // POSTCONDITION:
1528 : // - MinPoint%X <= XNext <= MaxPoint%X
1529 : // - LowerPoint%X < XNext < UpperPoint%X
1530 :
1531 : // Return value
1532 : bool BracketRoot;
1533 :
1534 : // Cannot use Secant method unless there are at least 2 points
1535 : // Also do not use Secant method more than once, i.e. NumHistory==3, in order to avoid
1536 : // the pathological case whereby the secant method always comes up short of bracketing
1537 : // the root because the function slope flattens as we come closer to either min/max point.
1538 5794018 : if (RootFinderData.NumHistory != 2) {
1539 4060206 : BracketRoot = false;
1540 4060206 : return BracketRoot;
1541 : }
1542 :
1543 : // Should not use Secant method if the last 2 points produced a warning
1544 1733812 : if (RootFinderData.StatusFlag == RootFinderStatus::WarningSingular || RootFinderData.StatusFlag == RootFinderStatus::WarningNonMonotonic) {
1545 5822 : BracketRoot = false;
1546 5822 : return BracketRoot;
1547 : }
1548 :
1549 : // Try to compute next root candidate using Secant formula
1550 1727990 : if (SecantFormula(RootFinderData, XNext)) {
1551 :
1552 : // Check that next candidate is consistent with min/max constraints and lower/upper brackets
1553 1727905 : if (CheckRootFinderCandidate(RootFinderData, XNext)) {
1554 1727905 : BracketRoot = true;
1555 1727905 : return BracketRoot;
1556 : }
1557 : }
1558 :
1559 85 : BracketRoot = false;
1560 :
1561 85 : return BracketRoot;
1562 : }
1563 :
1564 1025505 : Real64 BisectionMethod(RootFinderDataType &RootFinderData) // Data used by root finding algorithm
1565 : {
1566 : // FUNCTION INFORMATION:
1567 : // AUTHOR Dimitri Curtil (LBNL)
1568 : // DATE WRITTEN February 2006
1569 : // MODIFIED
1570 : // RE-ENGINEERED na
1571 :
1572 : // PURPOSE OF THIS FUNCTION:
1573 : // This function computes the next iterate using the bisection method (aka interval halving).
1574 : // Convergence rate is at best linear.
1575 : // PRECONDITION:
1576 : // Lower and upper points must be defined and distinct.
1577 : // POSTCONDITION:
1578 : // - LowerPoint%X < XCandidate < UpperPoint%X
1579 : // - RootFinderData%CurrentMethodType update with current solution method.
1580 :
1581 : // Return value
1582 : Real64 BisectionMethod;
1583 :
1584 1025505 : RootFinderData.CurrentMethodType = RootFinderMethod::Bisection;
1585 1025505 : BisectionMethod = (RootFinderData.LowerPoint.X + RootFinderData.UpperPoint.X) / 2.0;
1586 :
1587 1025505 : return BisectionMethod;
1588 : }
1589 :
1590 56421 : Real64 FalsePositionMethod(RootFinderDataType &RootFinderData) // Data used by root finding algorithm
1591 : {
1592 : // FUNCTION INFORMATION:
1593 : // AUTHOR Dimitri Curtil (LBNL)
1594 : // DATE WRITTEN February 2006
1595 : // MODIFIED
1596 : // RE-ENGINEERED na
1597 :
1598 : // PURPOSE OF THIS FUNCTION:
1599 : // This function computes the next iterate using the false position method (aka regula falsi).
1600 : // If new iterate does not lie within the lower and upper points then
1601 : // the Bisection method is used instead.
1602 : // Convergence rate is at best superlinear.
1603 : // PRECONDITION:
1604 : // Lower and upper points must be defined and distinct.
1605 : // POSTCONDITION:
1606 : // - LowerPoint%X < XCandidate < UpperPoint%X
1607 : // - RootFinderData%CurrentMethodType update with current solution method.
1608 :
1609 : // Return value
1610 : Real64 FalsePositionMethod;
1611 :
1612 : // FUNCTION LOCAL VARIABLE DECLARATIONS:
1613 : Real64 XCandidate;
1614 : Real64 Num;
1615 : Real64 Den;
1616 :
1617 56421 : Num = RootFinderData.UpperPoint.X - RootFinderData.LowerPoint.X;
1618 56421 : Den = RootFinderData.UpperPoint.Y - RootFinderData.LowerPoint.Y;
1619 :
1620 56421 : if (Den != 0.0) {
1621 : // False position method
1622 56421 : RootFinderData.CurrentMethodType = RootFinderMethod::FalsePosition;
1623 56421 : XCandidate = RootFinderData.LowerPoint.X - RootFinderData.LowerPoint.Y * Num / Den;
1624 :
1625 : // Check that new candidate is within range and brackets
1626 56421 : if (!CheckRootFinderCandidate(RootFinderData, XCandidate)) {
1627 : // Recovery method
1628 0 : XCandidate = BisectionMethod(RootFinderData);
1629 : }
1630 : } else {
1631 : // Recovery method
1632 0 : XCandidate = BisectionMethod(RootFinderData);
1633 : }
1634 :
1635 56421 : FalsePositionMethod = XCandidate;
1636 56421 : return FalsePositionMethod;
1637 : }
1638 :
1639 1371909 : Real64 SecantMethod(RootFinderDataType &RootFinderData) // Data used by root finding algorithm
1640 : {
1641 : // FUNCTION INFORMATION:
1642 : // AUTHOR Dimitri Curtil (LBNL)
1643 : // DATE WRITTEN February 2006
1644 : // MODIFIED
1645 : // RE-ENGINEERED na
1646 :
1647 : // PURPOSE OF THIS FUNCTION:
1648 : // This function computes the next iterate using the secant method.
1649 : // If new iterate does not lie within the lower and upper points then
1650 : // the false position method is used instead.
1651 : // Convergence rate is at best superlinear.
1652 : // PRECONDITION:
1653 : // There must be at least 2 history points so that RootFinderData%Increment is defined.
1654 : // See FUNCTION SecantFormula.
1655 : // POSTCONDITION:
1656 : // - LowerPoint%X < XCandidate < UpperPoint%X
1657 : // - RootFinderData%CurrentMethodType update with current solution method.
1658 :
1659 : // Return value
1660 : Real64 SecantMethod;
1661 :
1662 : // FUNCTION LOCAL VARIABLE DECLARATIONS:
1663 : Real64 XCandidate;
1664 :
1665 : // Recover with false position
1666 1371909 : if (SecantFormula(RootFinderData, XCandidate)) {
1667 : // Secant method
1668 1371909 : RootFinderData.CurrentMethodType = RootFinderMethod::Secant;
1669 :
1670 : // Check that new candidate is within range and brackets
1671 1371909 : if (!CheckRootFinderCandidate(RootFinderData, XCandidate)) {
1672 : // Recovery method
1673 0 : XCandidate = FalsePositionMethod(RootFinderData);
1674 : }
1675 : } else {
1676 : // Recovery method
1677 0 : XCandidate = FalsePositionMethod(RootFinderData);
1678 : }
1679 :
1680 1371909 : SecantMethod = XCandidate;
1681 1371909 : return SecantMethod;
1682 : }
1683 :
1684 3099899 : bool SecantFormula(RootFinderDataType const &RootFinderData, // Data used by root finding algorithm
1685 : Real64 &XNext // Result from Secant formula if possible to compute
1686 : )
1687 : {
1688 : // FUNCTION INFORMATION:
1689 : // AUTHOR Dimitri Curtil (LBNL)
1690 : // DATE WRITTEN April 2006
1691 : // MODIFIED
1692 : // RE-ENGINEERED na
1693 :
1694 : // PURPOSE OF THIS FUNCTION:
1695 : // This function computes the next iterate using the secant formula.
1696 : // If the new iterate cannot be computed the function returns FALSE, else TRUE.
1697 : // Convergence rate is at best superlinear.
1698 : // PRECONDITION:
1699 : // There must be at least 2 history points so that RootFinderData%Increment is defined.
1700 : // POSTCONDITION:
1701 : // XNext contains the result from applying the Secant formula.
1702 : // If XNext could not be computed then leave XNext unchanged.
1703 :
1704 : // Return value
1705 : bool SecantFormula;
1706 :
1707 : // FUNCTION LOCAL VARIABLE DECLARATIONS:
1708 : Real64 Num;
1709 : Real64 Den;
1710 :
1711 3099899 : Num = RootFinderData.Increment.X;
1712 3099899 : Den = RootFinderData.Increment.Y;
1713 :
1714 : // Cannot use secant with infinite slope (Den==0).
1715 : // Cannot use secant with null slope (Num==0).
1716 3099899 : if (Den != 0.0 && Num != 0.0) {
1717 3099814 : XNext = RootFinderData.CurrentPoint.X - RootFinderData.CurrentPoint.Y * Num / Den;
1718 3099814 : SecantFormula = true;
1719 : } else {
1720 85 : SecantFormula = false;
1721 : }
1722 :
1723 3099899 : return SecantFormula;
1724 : }
1725 :
1726 4976087 : Real64 BrentMethod(RootFinderDataType &RootFinderData) // Data used by root finding algorithm
1727 : {
1728 : // FUNCTION INFORMATION:
1729 : // AUTHOR Dimitri Curtil (LBNL)
1730 : // DATE WRITTEN March 2006
1731 : // MODIFIED
1732 : // RE-ENGINEERED na
1733 :
1734 : // PURPOSE OF THIS FUNCTION:
1735 : // This function computes the next iterate using the Brent's method.
1736 : // If new iterate does not lie within the lower and upper points then
1737 : // the secant method is used instead.
1738 : // Convergence rate is at best quadratic.
1739 : // PRECONDITION:
1740 : // Lower and upper points must be defined and distinct.
1741 : // POSTCONDITION:
1742 : // - LowerPoint%X < XCandidate < UpperPoint%X
1743 : // - RootFinderData%CurrentMethodType update with current solution method.
1744 :
1745 : // METHODOLOGY EMPLOYED:
1746 : // Inverse quadratic interpolation using the last 3 best iterates.
1747 : // The next root estimate is x = B + P/Q whereby B is the current best estimate
1748 : // of the root.
1749 :
1750 : // Return value
1751 : Real64 BrentMethod;
1752 :
1753 : // FUNCTION LOCAL VARIABLE DECLARATIONS:
1754 : Real64 XCandidate;
1755 : Real64 A;
1756 : Real64 FA;
1757 : Real64 B;
1758 : Real64 FB;
1759 : Real64 C;
1760 : Real64 FC;
1761 : Real64 R;
1762 : Real64 S;
1763 : Real64 T;
1764 : Real64 P;
1765 : Real64 Q;
1766 :
1767 : // Only attempt Brent's method if enough history points are available
1768 : // and if the root finder is converging (not diverging) since the previous
1769 : // iterate.
1770 : // We assume that;
1771 : // - the root is bracketed between the lower and upper points (see AdvanceRootFinder() ).
1772 : // - there are at least 3 history points
1773 4976087 : if (RootFinderData.NumHistory == 3) {
1774 :
1775 3604178 : A = RootFinderData.History(2).X;
1776 3604178 : FA = RootFinderData.History(2).Y;
1777 3604178 : B = RootFinderData.History(1).X;
1778 3604178 : FB = RootFinderData.History(1).Y;
1779 3604178 : C = RootFinderData.History(3).X;
1780 3604178 : FC = RootFinderData.History(3).Y;
1781 :
1782 : // Should never happen if CheckRootFinderConvergence() is invoked prior to this subroutine
1783 3604178 : if (FC == 0.0) {
1784 0 : BrentMethod = C;
1785 0 : return BrentMethod;
1786 : // Should never happen if CheckRootFinderConvergence() is invoked prior to this subroutine
1787 3604178 : } else if (FA == 0.0) {
1788 0 : BrentMethod = A;
1789 0 : return BrentMethod;
1790 : } else {
1791 3604178 : R = FB / FC;
1792 3604178 : S = FB / FA;
1793 3604178 : T = FA / FC;
1794 :
1795 3604178 : P = S * (T * (R - T) * (C - B) - (1.0 - R) * (B - A));
1796 3604178 : Q = (T - 1.0) * (R - 1.0) * (S - 1.0);
1797 :
1798 : // Only accept correction if it is small enough (75% of previous increment)
1799 3604178 : if (std::abs(P) <= 0.75 * std::abs(Q * RootFinderData.Increment.X)) {
1800 2578673 : RootFinderData.CurrentMethodType = RootFinderMethod::Brent;
1801 2578673 : XCandidate = B + P / Q;
1802 :
1803 : // Check that new candidate is within range and brackets
1804 2578673 : if (!CheckRootFinderCandidate(RootFinderData, XCandidate)) {
1805 : // Recovery method
1806 3436 : XCandidate = FalsePositionMethod(RootFinderData);
1807 : }
1808 : } else {
1809 : // Recover from bad correction with bisection
1810 : // Bisection produced the best numerical performance in testing compared to
1811 : // - Secant
1812 : // - False position (very slow recovery)
1813 1025505 : XCandidate = BisectionMethod(RootFinderData);
1814 : }
1815 : }
1816 : } else {
1817 : // Not enough history to try Brent's method yet: use Secant's method
1818 1371909 : XCandidate = SecantMethod(RootFinderData);
1819 : }
1820 :
1821 4976087 : BrentMethod = XCandidate;
1822 4976087 : return BrentMethod;
1823 : }
1824 :
1825 0 : void WriteRootFinderTraceHeader(InputOutputFile &TraceFile) // Unit for trace file
1826 : {
1827 : // SUBROUTINE INFORMATION:
1828 : // AUTHOR Dimitri Curtil (LBNL)
1829 : // DATE WRITTEN March 2006
1830 : // MODIFIED
1831 : // RE-ENGINEERED na
1832 :
1833 : // PURPOSE OF THIS SUBROUTINE:
1834 : // This subroutine writes the header for the trace file to the specified
1835 : // file unit using CSV formatting.
1836 :
1837 0 : print(TraceFile,
1838 : "{},{},{},{},{},{},{},{},{},{},{},{},{},{},{},{},{},{},{},{},",
1839 : "Status",
1840 : "Method",
1841 : "CurrentPoint%X",
1842 : "CurrentPoint%Y",
1843 : "XCandidate",
1844 : "ConvergenceRate",
1845 : "MinPoint%X",
1846 : "MinPoint%Y",
1847 : "LowerPoint%X",
1848 : "LowerPoint%Y",
1849 : "UpperPoint%X",
1850 : "UpperPoint%Y",
1851 : "MaxPoint%X",
1852 : "MaxPoint%Y",
1853 : "History(1)%X",
1854 : "History(1)%Y",
1855 : "History(2)%X",
1856 : "History(2)%Y",
1857 : "History(3)%X",
1858 : "History(3)%Y");
1859 0 : }
1860 :
1861 0 : void WriteRootFinderTrace(InputOutputFile &TraceFile, // Unit for trace file
1862 : RootFinderDataType const &RootFinderData // Data used by root finding algorithm
1863 : )
1864 : {
1865 : // SUBROUTINE INFORMATION:
1866 : // AUTHOR Dimitri Curtil (LBNL)
1867 : // DATE WRITTEN March 2006
1868 : // MODIFIED
1869 : // RE-ENGINEERED na
1870 :
1871 : // PURPOSE OF THIS SUBROUTINE:
1872 : // This subroutine writes the current state of the root finder data to the trace file
1873 : // unit using CSV formatting.
1874 :
1875 0 : print(TraceFile, "{},{},", RootFinderData.StatusFlag, RootFinderData.CurrentMethodType);
1876 :
1877 : // Only show current point if defined.
1878 0 : WritePoint(TraceFile, RootFinderData.CurrentPoint, false);
1879 :
1880 0 : print(TraceFile, "{:20.10F},{:20.10F},", RootFinderData.XCandidate, RootFinderData.ConvergenceRate);
1881 :
1882 : // Always show min and max points.
1883 : // Only show lower and upper points if defined.
1884 0 : WritePoint(TraceFile, RootFinderData.MinPoint, true);
1885 0 : WritePoint(TraceFile, RootFinderData.LowerPoint, false);
1886 0 : WritePoint(TraceFile, RootFinderData.UpperPoint, false);
1887 0 : WritePoint(TraceFile, RootFinderData.MaxPoint, true);
1888 : // Only show history points if defined.
1889 0 : WritePoint(TraceFile, RootFinderData.History(1), false);
1890 0 : WritePoint(TraceFile, RootFinderData.History(2), false);
1891 0 : WritePoint(TraceFile, RootFinderData.History(3), false);
1892 0 : }
1893 :
1894 0 : void WritePoint(InputOutputFile &TraceFile, // Unit for trace file
1895 : PointType const &PointData, // Point data structure
1896 : bool const ShowXValue)
1897 : {
1898 : // SUBROUTINE INFORMATION:
1899 : // AUTHOR Dimitri Curtil (LBNL)
1900 : // DATE WRITTEN March 2006
1901 : // MODIFIED
1902 : // RE-ENGINEERED na
1903 :
1904 : // PURPOSE OF THIS SUBROUTINE:
1905 : // This subroutine writes the current point data to the trace file
1906 : // unit using CSV formatting.
1907 : // If not defined writes an empty string instead.
1908 :
1909 0 : if (PointData.DefinedFlag) {
1910 0 : print(TraceFile, "{:20.10F},{:20.10F},", PointData.X, PointData.Y);
1911 : } else {
1912 0 : if (ShowXValue) {
1913 0 : print(TraceFile, "{:20.10F},,", PointData.X);
1914 : } else {
1915 0 : print(TraceFile, ",,");
1916 : }
1917 : }
1918 0 : }
1919 :
1920 0 : void DebugRootFinder(InputOutputFile &DebugFile, // File unit where to write debugging info
1921 : RootFinderDataType const &RootFinderData // Data used by root finding algorithm
1922 : )
1923 : {
1924 : // SUBROUTINE INFORMATION:
1925 : // AUTHOR Dimitri Curtil (LBNL)
1926 : // DATE WRITTEN April 2006
1927 : // MODIFIED
1928 : // RE-ENGINEERED na
1929 :
1930 : // PURPOSE OF THIS SUBROUTINE:
1931 : // This subroutine writes the current min/max range and lower/upper bracket to
1932 : // the standard output file.
1933 : // Used only for debugging.
1934 :
1935 0 : print(DebugFile, "Current = ");
1936 0 : WritePoint(DebugFile, RootFinderData.CurrentPoint, true);
1937 0 : print(DebugFile, "\n");
1938 :
1939 0 : print(DebugFile, "Min = ");
1940 0 : WritePoint(DebugFile, RootFinderData.MinPoint, true);
1941 0 : print(DebugFile, "\n");
1942 :
1943 0 : print(DebugFile, "Lower = ");
1944 0 : WritePoint(DebugFile, RootFinderData.LowerPoint, false);
1945 0 : print(DebugFile, "\n");
1946 :
1947 0 : print(DebugFile, "Upper = ");
1948 0 : WritePoint(DebugFile, RootFinderData.UpperPoint, false);
1949 0 : print(DebugFile, "\n");
1950 :
1951 0 : print(DebugFile, "Max = ");
1952 0 : WritePoint(DebugFile, RootFinderData.MaxPoint, true);
1953 0 : print(DebugFile, "\n");
1954 0 : }
1955 :
1956 0 : void WriteRootFinderStatus(InputOutputFile &File, // File unit where to write the status description
1957 : RootFinderDataType const &RootFinderData // Data used by root finding algorithm
1958 : )
1959 : {
1960 : // SUBROUTINE INFORMATION:
1961 : // AUTHOR Dimitri Curtil (LBNL)
1962 : // DATE WRITTEN May 2006
1963 : // MODIFIED
1964 : // RE-ENGINEERED na
1965 :
1966 0 : switch (RootFinderData.StatusFlag) {
1967 0 : case RootFinderStatus::OK: {
1968 0 : print(File, "Found unconstrained root");
1969 0 : } break;
1970 0 : case RootFinderStatus::OKMin: {
1971 0 : print(File, "Found min constrained root");
1972 0 : } break;
1973 0 : case RootFinderStatus::OKMax: {
1974 0 : print(File, "Found max constrained root");
1975 0 : } break;
1976 0 : case RootFinderStatus::OKRoundOff: {
1977 0 : print(File, "Detected round-off convergence in bracket");
1978 0 : } break;
1979 0 : case RootFinderStatus::WarningSingular: {
1980 0 : print(File, "Detected singularity warning");
1981 0 : } break;
1982 0 : case RootFinderStatus::WarningNonMonotonic: {
1983 0 : print(File, "Detected non-monotonicity warning");
1984 0 : } break;
1985 0 : case RootFinderStatus::ErrorRange: {
1986 0 : print(File, "Detected out-of-range error");
1987 0 : } break;
1988 0 : case RootFinderStatus::ErrorBracket: {
1989 0 : print(File, "Detected bracket error");
1990 0 : } break;
1991 0 : case RootFinderStatus::ErrorSlope: {
1992 0 : print(File, "Detected slope error");
1993 0 : } break;
1994 0 : case RootFinderStatus::ErrorSingular: {
1995 0 : print(File, "Detected singularity error");
1996 0 : } break;
1997 0 : default: {
1998 0 : print(File, "Detected bad root finder status");
1999 0 : } break;
2000 : }
2001 0 : }
2002 :
2003 : } // namespace EnergyPlus::RootFinder
|
