1#  ___________________________________________________________________________
2#
3#  Pyomo: Python Optimization Modeling Objects
4#  Copyright 2017 National Technology and Engineering Solutions of Sandia, LLC
5#  Under the terms of Contract DE-NA0003525 with National Technology and
6#  Engineering Solutions of Sandia, LLC, the U.S. Government retains certain
7#  rights in this software.
8#  This software is distributed under the 3-clause BSD License.
9#  ___________________________________________________________________________
10
11import math
12import logging
13from pyomo.common.errors import DeveloperError, InfeasibleConstraintException, PyomoException
14
15logger = logging.getLogger(__name__)
16inf = float('inf')
17
18
19class IntervalException(PyomoException):
20    pass
21
22
23def add(xl, xu, yl, yu):
24    return xl + yl, xu + yu
25
26
27def sub(xl, xu, yl, yu):
28    return xl - yu, xu - yl
29
30
31def mul(xl, xu, yl, yu):
32    options = [xl*yl, xl*yu, xu*yl, xu*yu]
33    if any(math.isnan(i) for i in options):
34        lb = -inf
35        ub = inf
36    else:
37        lb = min(options)
38        ub = max(options)
39    return lb, ub
40
41
42def inv(xl, xu, feasibility_tol):
43    """
44    The case where xl is very slightly positive but should be very slightly negative (or xu is very slightly negative
45    but should be very slightly positive) should not be an issue. Suppose xu is 2 and xl is 1e-15 but should be -1e-15.
46    The bounds obtained from this function will be [0.5, 1e15] or [0.5, inf), depending on the value of
47    feasibility_tol. The true bounds are (-inf, -1e15] U [0.5, inf), where U is union. The exclusion of (-inf, -1e15]
48    should be acceptable. Additionally, it very important to return a non-negative interval when xl is non-negative.
49    """
50    if xu - xl <= -feasibility_tol:
51        raise InfeasibleConstraintException(f'lower bound is greater than upper bound in inv; xl: {xl}; xu: {xu}')
52    elif xu <= 0 <= xl:
53        # This has to return -inf to inf because it could later be multiplied by 0
54        lb = -inf
55        ub = inf
56    elif xl < 0 < xu:
57        lb = -inf
58        ub = inf
59    elif 0 <= xl <= feasibility_tol:
60        # xu must be strictly positive
61        ub = inf
62        lb = 1.0 / xu
63    elif xl > feasibility_tol:
64        # xl and xu must be strictly positive
65        ub = 1.0 / xl
66        lb = 1.0 / xu
67    elif -feasibility_tol <= xu <= 0:
68        # xl must be strictly negative
69        lb = -inf
70        ub = 1.0 / xl
71    elif xu < -feasibility_tol:
72        # xl and xu must be strictly negative
73        ub = 1.0 / xl
74        lb = 1.0 / xu
75    else:
76        # everything else
77        lb = -inf
78        ub = inf
79    return lb, ub
80
81
82def div(xl, xu, yl, yu, feasibility_tol):
83    lb, ub = mul(xl, xu, *inv(yl, yu, feasibility_tol))
84    return lb, ub
85
86
87def power(xl, xu, yl, yu, feasibility_tol):
88    """
89    Compute bounds on x**y.
90    """
91    if xl > 0:
92        """
93        If x is always positive, things are simple. We only need to worry about the sign of y.
94        """
95        if yl < 0 < yu:
96            lb = min(xu ** yl, xl ** yu)
97            ub = max(xl ** yl, xu ** yu)
98        elif yl >= 0:
99            lb = min(xl**yl, xl**yu)
100            ub = max(xu**yl, xu**yu)
101        else:  # yu <= 0:
102            lb = min(xu**yl, xu**yu)
103            ub = max(xl**yl, xl**yu)
104    elif xl == 0:
105        if yl >= 0:
106            lb = min(xl ** yl, xl ** yu)
107            ub = max(xu ** yl, xu ** yu)
108        elif yu <= 0:
109            lb, ub = inv(*power(xl, xu, *sub(0, 0, yl, yu), feasibility_tol), feasibility_tol)
110        else:
111            lb1, ub1 = power(xl, xu, 0, yu, feasibility_tol)
112            lb2, ub2 = power(xl, xu, yl, 0, feasibility_tol)
113            lb = min(lb1, lb2)
114            ub = max(ub1, ub2)
115    elif yl == yu and yl == round(yl):
116        # the exponent is an integer, so x can be negative
117        """
118        The logic here depends on several things:
119        1) The sign of x
120        2) The sign of y
121        3) Whether y is even or odd.
122
123        There are also special cases to avoid math domain errors.
124        """
125        y = yl
126        if xu <= 0:
127            if y < 0:
128                if y % 2 == 0:
129                    lb = xl ** y
130                    if xu == 0:
131                        ub = inf
132                    else:
133                        ub = xu ** y
134                else:
135                    if xu == 0:
136                        lb = -inf
137                        ub = inf
138                    else:
139                        lb = xu ** y
140                        ub = xl ** y
141            else:
142                if y % 2 == 0:
143                    lb = xu ** y
144                    ub = xl ** y
145                else:
146                    lb = xl ** y
147                    ub = xu ** y
148        else:
149            if y < 0:
150                if y % 2 == 0:
151                    lb = min(xl ** y, xu ** y)
152                    ub = inf
153                else:
154                    lb = - inf
155                    ub = inf
156            else:
157                if y % 2 == 0:
158                    lb = 0
159                    ub = max(xl ** y, xu ** y)
160                else:
161                    lb = xl ** y
162                    ub = xu ** y
163    elif yl == yu:
164        # the exponent has to be fractional, so x must be positive
165        if xu < 0:
166            msg = 'Cannot raise a negative number to the power of {0}.\n'.format(yl)
167            msg += 'The upper bound of a variable raised to the power of {0} is {1}'.format(yl, xu)
168            raise InfeasibleConstraintException(msg)
169        xl = 0
170        lb, ub = power(xl, xu, yl, yu, feasibility_tol)
171    else:
172        lb = -inf
173        ub = inf
174
175    return lb, ub
176
177
178def _inverse_power1(zl, zu, yl, yu, orig_xl, orig_xu, feasibility_tol):
179    """
180    z = x**y => compute bounds on x.
181
182    First, start by computing bounds on x with
183
184        x = exp(ln(z) / y)
185
186    However, if y is an integer, then x can be negative, so there are several special cases. See the docs below.
187    """
188    xl, xu = log(zl, zu)
189    xl, xu = div(xl, xu, yl, yu, feasibility_tol)
190    xl, xu = exp(xl, xu)
191
192    # if y is an integer, then x can be negative
193    if yl == yu and yl == round(yl):  # y is a fixed integer
194        y = yl
195        if y == 0:
196            # Anything to the power of 0 is 1, so if y is 0, then x can be anything
197            # (assuming zl <= 1 <= zu, which is enforced when traversing the tree in the other direction)
198            xl = -inf
199            xu = inf
200        elif y % 2 == 0:
201            """
202            if y is even, then there are two primary cases (note that it is much easier to walk through these
203            while looking at plots):
204            case 1: y is positive
205                x**y is convex, positive, and symmetric. The bounds on x depend on the lower bound of z. If zl <= 0,
206                then xl should simply be -xu. However, if zl > 0, then we may be able to say something better. For
207                example, if the original lower bound on x is positive, then we can keep xl computed from
208                x = exp(ln(z) / y). Furthermore, if the original lower bound on x is larger than -xl computed from
209                x = exp(ln(z) / y), then we can still keep the xl computed from x = exp(ln(z) / y). Similar logic
210                applies to the upper bound of x.
211            case 2: y is negative
212                The ideas are similar to case 1.
213            """
214            if zu + feasibility_tol < 0:
215                raise InfeasibleConstraintException('Infeasible. Anything to the power of {0} must be positive.'.format(y))
216            if y > 0:
217                if zu <= 0:
218                    _xl = 0
219                    _xu = 0
220                elif zl <= 0:
221                    _xl = -xu
222                    _xu = xu
223                else:
224                    if orig_xl <= -xl + feasibility_tol:
225                        _xl = -xu
226                    else:
227                        _xl = xl
228                    if orig_xu < xl - feasibility_tol:
229                        _xu = -xl
230                    else:
231                        _xu = xu
232                xl = _xl
233                xu = _xu
234            else:
235                if zu == 0:
236                    raise InfeasibleConstraintException('Infeasible. Anything to the power of {0} must be positive.'.format(y))
237                elif zl <= 0:
238                    _xl = -inf
239                    _xu = inf
240                else:
241                    if orig_xl <= -xl + feasibility_tol:
242                        _xl = -xu
243                    else:
244                        _xl = xl
245                    if orig_xu < xl - feasibility_tol:
246                        _xu = -xl
247                    else:
248                        _xu = xu
249                xl = _xl
250                xu = _xu
251        else:  # y % 2 == 1
252            """
253            y is odd.
254            Case 1: y is positive
255                x**y is monotonically increasing. If y is positive, then we can can compute the bounds on x using
256                x = z**(1/y) and the signs on xl and xu depend on the signs of zl and zu.
257            Case 2: y is negative
258                Again, this is easier to visualize with a plot. x**y approaches zero when x approaches -inf or inf.
259                Thus, if zl < 0 < zu, then no bounds can be inferred for x. If z is positive (zl >=0 ) then we can
260                use the bounds computed from x = exp(ln(z) / y). If z is negative (zu <= 0), then we live in the
261                bottom left quadrant, xl depends on zu, and xu depends on zl.
262            """
263            if y > 0:
264                xl = abs(zl)**(1.0/y)
265                xl = math.copysign(xl, zl)
266                xu = abs(zu)**(1.0/y)
267                xu = math.copysign(xu, zu)
268            else:
269                if zl >= 0:
270                    pass
271                elif zu <= 0:
272                    if zu == 0:
273                        xl = -inf
274                    else:
275                        xl = -abs(zu)**(1.0/y)
276                    if zl == 0:
277                        xu = -inf
278                    else:
279                        xu = -abs(zl)**(1.0/y)
280                else:
281                    xl = -inf
282                    xu = inf
283
284    return xl, xu
285
286
287def _inverse_power2(zl, zu, xl, xu, feasiblity_tol):
288    """
289    z = x**y => compute bounds on y
290    y = ln(z) / ln(x)
291
292    This function assumes the exponent can be fractional, so x must be positive. This method should not be called
293    if the exponent is an integer.
294    """
295    if xu <= 0:
296        raise IntervalException('Cannot raise a negative variable to a fractional power.')
297    if (xl > 0 and zu <= 0) or (xl >= 0 and zu < 0):
298        raise InfeasibleConstraintException('A positive variable raised to the power of anything must be positive.')
299    lba, uba = log(zl, zu)
300    lbb, ubb = log(xl, xu)
301    yl, yu = div(lba, uba, lbb, ubb, feasiblity_tol)
302    return yl, yu
303
304
305def exp(xl, xu):
306    try:
307        lb = math.exp(xl)
308    except OverflowError:
309        lb = math.inf
310    try:
311        ub = math.exp(xu)
312    except OverflowError:
313        ub = math.inf
314    return lb, ub
315
316
317def log(xl, xu):
318    if xl > 0:
319        lb = math.log(xl)
320    else:
321        lb = -inf
322    if xu > 0:
323        ub = math.log(xu)
324    else:
325        ub = -inf
326    return lb, ub
327
328
329def log10(xl, xu):
330    if xl > 0:
331        lb = math.log10(xl)
332    else:
333        lb = -inf
334    if xu > 0:
335        ub = math.log10(xu)
336    else:
337        ub = -inf
338    return lb, ub
339
340
341def sin(xl, xu):
342    """
343
344    Parameters
345    ----------
346    xl: float
347    xu: float
348
349    Returns
350    -------
351    lb: float
352    ub: float
353    """
354
355    # if there is a minimum between xl and xu, then the lower bound is -1. Minimums occur at 2*pi*n - pi/2
356    # find the minimum value of i such that 2*pi*i - pi/2 >= xl. Then round i up. If 2*pi*i - pi/2 is still less
357    # than or equal to xu, then there is a minimum between xl and xu. Thus the lb is -1. Otherwise, the minimum
358    # occurs at either xl or xu
359    if xl <= -inf or xu >= inf:
360        return -1, 1
361    pi = math.pi
362    i = (xl + pi / 2) / (2 * pi)
363    i = math.ceil(i)
364    x_at_min = 2 * pi * i - pi / 2
365    if x_at_min <= xu:
366        lb = -1
367    else:
368        lb = min(math.sin(xl), math.sin(xu))
369
370    # if there is a maximum between xl and xu, then the upper bound is 1. Maximums occur at 2*pi*n + pi/2
371    i = (xu - pi / 2) / (2 * pi)
372    i = math.floor(i)
373    x_at_max = 2 * pi * i + pi / 2
374    if x_at_max >= xl:
375        ub = 1
376    else:
377        ub = max(math.sin(xl), math.sin(xu))
378
379    return lb, ub
380
381
382def cos(xl, xu):
383    """
384
385    Parameters
386    ----------
387    xl: float
388    xu: float
389
390    Returns
391    -------
392    lb: float
393    ub: float
394    """
395
396    # if there is a minimum between xl and xu, then the lower bound is -1. Minimums occur at 2*pi*n - pi
397    # find the minimum value of i such that 2*pi*i - pi >= xl. Then round i up. If 2*pi*i - pi/2 is still less
398    # than or equal to xu, then there is a minimum between xl and xu. Thus the lb is -1. Otherwise, the minimum
399    # occurs at either xl or xu
400    if xl <= -inf or xu >= inf:
401        return -1, 1
402    pi = math.pi
403    i = (xl + pi) / (2 * pi)
404    i = math.ceil(i)
405    x_at_min = 2 * pi * i - pi
406    if x_at_min <= xu:
407        lb = -1
408    else:
409        lb = min(math.cos(xl), math.cos(xu))
410
411    # if there is a maximum between xl and xu, then the upper bound is 1. Maximums occur at 2*pi*n
412    i = (xu) / (2 * pi)
413    i = math.floor(i)
414    x_at_max = 2 * pi * i
415    if x_at_max >= xl:
416        ub = 1
417    else:
418        ub = max(math.cos(xl), math.cos(xu))
419
420    return lb, ub
421
422
423def tan(xl, xu):
424    """
425
426    Parameters
427    ----------
428    xl: float
429    xu: float
430
431    Returns
432    -------
433    lb: float
434    ub: float
435    """
436
437    # tan goes to -inf and inf at every pi*i + pi/2 (integer i). If one of these values is between xl and xu, then
438    # the lb is -inf and the ub is inf. Otherwise the minimum occurs at xl and the maximum occurs at xu.
439    # find the minimum value of i such that pi*i + pi/2 >= xl. Then round i up. If pi*i + pi/2 is still less
440    # than or equal to xu, then there is an undefined point between xl and xu.
441    if xl <= -inf or xu >= inf:
442        return -inf, inf
443    pi = math.pi
444    i = (xl - pi / 2) / (pi)
445    i = math.ceil(i)
446    x_at_undefined = pi * i + pi / 2
447    if x_at_undefined <= xu:
448        lb = -inf
449        ub = inf
450    else:
451        lb = math.tan(xl)
452        ub = math.tan(xu)
453
454    return lb, ub
455
456
457def asin(xl, xu, yl, yu, feasibility_tol):
458    """
459    y = asin(x); propagate bounds from x to y
460    x = sin(y)
461
462    Parameters
463    ----------
464    xl: float
465    xu: float
466    yl: float
467    yu: float
468
469    Returns
470    -------
471    lb: float
472    ub: float
473    """
474    if xl < -1:
475        xl = -1
476    if xu > 1:
477        xu = 1
478
479    pi = math.pi
480
481    if yl <= -inf:
482        lb = yl
483    elif xl <= math.sin(yl) <= xu:
484        # if sin(yl) >= xl then yl satisfies the bounds on x, and the lower bound of y cannot be improved
485        lb = yl
486    elif math.sin(yl) < xl:
487        """
488        we can only push yl up from its current value to the next lowest value such that xl = sin(y). In other words,
489        we need to
490
491        min y
492        s.t.
493            xl = sin(y)
494            y >= yl
495
496        globally.
497        """
498        # first find the next minimum of x = sin(y). Minimums occur at y = 2*pi*n - pi/2 for integer n.
499        i = (yl + pi / 2) / (2 * pi)
500        i1 = math.floor(i)
501        i2 = math.ceil(i)
502        i1 = 2 * pi * i1 - pi / 2
503        i2 = 2 * pi * i2 - pi / 2
504        # now find the next value of y such that xl = sin(y). This can be computed by a distance from the minimum (i).
505        y_tmp = math.asin(xl)  # this will give me a value between -pi/2 and pi/2
506        dist = y_tmp - (-pi / 2)  # this is the distance between the minimum of the sin function and a value that
507        # satisfies xl = sin(y)
508        lb1 = i1 + dist
509        lb2 = i2 + dist
510        if lb1 >= yl - feasibility_tol:
511            lb = lb1
512        else:
513            lb = lb2
514    else:
515        # sin(yl) > xu
516        i = (yl - pi / 2) / (2 * pi)
517        i1 = math.floor(i)
518        i2 = math.ceil(i)
519        i1 = 2 * pi * i1 + pi / 2
520        i2 = 2 * pi * i2 + pi / 2
521        y_tmp = math.asin(xu)
522        dist = pi / 2 - y_tmp
523        lb1 = i1 + dist
524        lb2 = i2 + dist
525        if lb1 >= yl - feasibility_tol:
526            lb = lb1
527        else:
528            lb = lb2
529
530    # use the same logic for the maximum
531    if yu >= inf:
532        ub = yu
533    elif xl <= math.sin(yu) <= xu:
534        ub = yu
535    elif math.sin(yu) > xu:
536        i = (yu - pi / 2) / (2 * pi)
537        i1 = math.ceil(i)
538        i2 = math.floor(i)
539        i1 = 2 * pi * i1 + pi / 2
540        i2 = 2 * pi * i2 + pi / 2
541        y_tmp = math.asin(xu)
542        dist = pi / 2 - y_tmp
543        ub1 = i1 - dist
544        ub2 = i2 - dist
545        if ub1 <= yu + feasibility_tol:
546            ub = ub1
547        else:
548            ub = ub2
549    else:
550        # math.sin(yu) < xl
551        i = (yu + pi / 2) / (2 * pi)
552        i1 = math.ceil(i)
553        i2 = math.floor(i)
554        i1 = 2 * pi * i1 - pi / 2
555        i2 = 2 * pi * i2 - pi / 2
556        y_tmp = math.asin(xl)
557        dist = y_tmp - (-pi / 2)
558        ub1 = i1 - dist
559        ub2 = i2 - dist
560        if ub1 <= yu + feasibility_tol:
561            ub = ub1
562        else:
563            ub = ub2
564
565    return lb, ub
566
567
568def acos(xl, xu, yl, yu, feasibility_tol):
569    """
570    y = acos(x); propagate bounds from x to y
571    x = cos(y)
572
573    Parameters
574    ----------
575    xl: float
576    xu: float
577    yl: float
578    yu: float
579
580    Returns
581    -------
582    lb: float
583    ub: float
584    """
585    if xl < -1:
586        xl = -1
587    if xu > 1:
588        xu = 1
589
590    pi = math.pi
591
592    if yl <= -inf:
593        lb = yl
594    elif xl <= math.cos(yl) <= xu:
595        # if xl <= cos(yl) <= xu then yl satisfies the bounds on x, and the lower bound of y cannot be improved
596        lb = yl
597    elif math.cos(yl) < xl:
598        """
599        we can only push yl up from its current value to the next lowest value such that xl = cos(y). In other words,
600        we need to
601
602        min y
603        s.t.
604            xl = cos(y)
605            y >= yl
606
607        globally.
608        """
609        # first find the next minimum of x = cos(y). Minimums occur at y = 2*pi*n - pi for integer n.
610        i = (yl + pi) / (2 * pi)
611        i1 = math.floor(i)
612        i2 = math.ceil(i)
613        i1 = 2 * pi * i1 - pi
614        i2 = 2 * pi * i2 - pi
615        # now find the next value of y such that xl = cos(y). This can be computed by a distance from the minimum (i).
616        y_tmp = math.acos(xl)  # this will give me a value between 0 and pi
617        dist = pi - y_tmp  # this is the distance between the minimum of the sin function and a value that
618        # satisfies xl = sin(y)
619        lb1 = i1 + dist
620        lb2 = i2 + dist
621        if lb1 >= yl - feasibility_tol:
622            lb = lb1
623        else:
624            lb = lb2
625    else:
626        # cos(yl) > xu
627        # first find the next maximum of x = cos(y).
628        i = yl / (2 * pi)
629        i1 = math.floor(i)
630        i2 = math.ceil(i)
631        i1 = 2 * pi * i1
632        i2 = 2 * pi * i2
633        y_tmp = math.acos(xu)
634        dist = y_tmp
635        lb1 = i1 + dist
636        lb2 = i2 + dist
637        if lb1 >= yl - feasibility_tol:
638            lb = lb1
639        else:
640            lb = lb2
641
642    # use the same logic for the maximum
643    if yu >= inf:
644        ub = yu
645    elif xl <= math.cos(yu) <= xu:
646        ub = yu
647    elif math.cos(yu) > xu:
648        i = yu / (2 * pi)
649        i1 = math.ceil(i)
650        i2 = math.floor(i)
651        i1 = 2 * pi * i1
652        i2 = 2 * pi * i2
653        y_tmp = math.acos(xu)
654        dist = y_tmp
655        ub1 = i1 - dist
656        ub2 = i2 - dist
657        if ub1 <= yu + feasibility_tol:
658            ub = ub1
659        else:
660            ub = ub2
661    else:
662        # math.cos(yu) < xl
663        i = (yu + pi) / (2 * pi)
664        i1 = math.ceil(i)
665        i2 = math.floor(i)
666        i1 = 2 * pi * i1 - pi
667        i2 = 2 * pi * i2 - pi
668        y_tmp = math.acos(xl)
669        dist = pi - y_tmp
670        ub1 = i1 - dist
671        ub2 = i2 - dist
672        if ub1 <= yu + feasibility_tol:
673            ub = ub1
674        else:
675            ub = ub2
676
677    return lb, ub
678
679
680def atan(xl, xu, yl, yu):
681    """
682    y = atan(x); propagate bounds from x to y
683    x = tan(y)
684
685    Parameters
686    ----------
687    xl: float
688    xu: float
689    yl: float
690    yu: float
691
692    Returns
693    -------
694    lb: float
695    ub: float
696    """
697
698    pi = math.pi
699
700    # tan goes to -inf and inf at every pi*i + pi/2 (integer i).
701    if xl <= -inf or yl <= -inf:
702        lb = yl
703    else:
704        i = (yl - pi / 2) / pi
705        i = math.floor(i)
706        i = pi * i + pi / 2
707        y_tmp = math.atan(xl)
708        dist = y_tmp - (-pi/2)
709        lb = i + dist
710
711    if xu >= inf or yu >= inf:
712        ub = yu
713    else:
714        i = (yu - pi / 2) / pi
715        i = math.ceil(i)
716        i = pi * i + pi / 2
717        y_tmp = math.atan(xu)
718        dist = pi / 2 - y_tmp
719        ub = i - dist
720
721    if yl > lb:
722        lb = yl
723    if yu < ub:
724        ub = yu
725
726    return lb, ub
727