User guide · 2.0 development

Advanced Features¶

Monomial and Signomial Terms¶

For geometric programming, you can use signomial terms:

problem = SHOTpy.Problem(solver)

x = SHOTpy.Variable("x", SHOTpy.VariableType.Real, 0.1, 10.0)
y = SHOTpy.Variable("y", SHOTpy.VariableType.Real, 0.1, 10.0)
problem.addVariable(x)
problem.addVariable(y)

# Create a signomial element: x^2 * y^(-1)
element1 = SHOTpy.SignomialElement(x, 2.0)
element2 = SHOTpy.SignomialElement(y, -1.0)

# Create signomial term: 3.0 * x^2 * y^(-1) - coefficient first, then elements
signomial = SHOTpy.SignomialTerm(3.0, [element1, element2])

# Alternative: Create monomial using variable-power pairs directly
monomial = SHOTpy.SignomialTerm(2.5, [(x, 2.0), (y, 3.0)])

# Use in constraint
constraint = SHOTpy.NonlinearConstraint("signomial_constraint", -SHOTpy.SHOT_DBL_MAX, 100.0)
constraint.add(signomial)
problem.addConstraint(constraint)

# Simple objective
objective = SHOTpy.LinearObjectiveFunction(SHOTpy.ObjectiveDirection.Minimize)
objective.add(SHOTpy.LinearTerm(1.0, x))
problem.setObjective(objective)

problem.finalize()
print("Problem with signomial terms created")
print(problem.toString())

Examining Solution Details¶

# Solve a simple problem
solver = SHOTpy.Solver()
problem = SHOTpy.Problem(solver)

x = problem.addVariable("x", SHOTpy.VariableType.Real, 0.0, 10.0)
y = problem.addVariable("y", SHOTpy.VariableType.Real, 0.0, 10.0)

problem.setObjective(x**2 + y**2, SHOTpy.ObjectiveDirection.Minimize)
problem.addConstraint(x + y >= 3, "sum")

problem.finalize()
assert solver.setProblem(problem)
assert solver.solveProblem()
assert solver.getModelReturnStatus() == SHOTpy.ModelReturnStatus.OptimalGlobal

# Get detailed results
print("Solution Details:")
sol = solver.getPrimalSolution()
print(f"  Primal objective: {sol.objValue:.6f}")
print(f"  Dual bound: {solver.getGlobalDualBound():.6f}")

solutions = solver.getPrimalSolutions()
print(f"  Number of solutions found: {len(solutions)}")

if len(solutions) > 0:
    print(f"  Best solution: x = {solutions[0].point[0]:.4f}, y = {solutions[0].point[1]:.4f}")

# Get problem name (if set)
print(f"  Problem name: {problem.name if hasattr(problem, 'name') else 'N/A'}")

# The optimum is x = y = 1.5, and the dual bound proves it up to the termination tolerances, an absolute or
# relative objective gap of 0.001
assert abs(sol.objValue - 4.5) < 1e-3
assert solver.getAbsoluteObjectiveGap() <= 1e-3 or solver.getRelativeObjectiveGap() <= 1e-3

Checking Problem Properties¶

# Create a problem with various features
problem = SHOTpy.Problem(solver)

x = problem.addVariable("x", SHOTpy.VariableType.Real, 0.1, 10.0)
y = problem.addVariable("y", SHOTpy.VariableType.Integer, 0.0, 5.0)

# Mix of linear, quadratic, and nonlinear terms
problem.setObjective(2 * x + y**2 + SHOTpy.log(x), SHOTpy.ObjectiveDirection.Minimize)

problem.addConstraint(x >= 1, "linear")
problem.addConstraint(SHOTpy.exp(x) <= 10, "nonlinear")

problem.finalize()

# Display problem structure with indicators
print("Problem Structure:")
print(problem.toString())
print("\nIndicators legend:")
print("  L = Linear terms")
print("  Q = Quadratic terms")
print("  S = Signomial terms (power functions)")
print("  M = Monomial terms")
print("  E = General nonlinear expressions")

Evaluating Functions¶

The objective function and the constraints can be evaluated at a point, which is a list with a value for each variable in the order of their indexes:

Method	Returns
`objectiveFunction.calculateValue(point)`	the value of the objective function, including its constant
`constraint.calculateFunctionValue(point)`	f(x) of the constraint L ≤ f(x) ≤ U
`constraint.calculateNumericValue(point)`	a `NumericConstraintValue` with `functionValue`, `isFulfilled`, `error` and the deviations from L and U
`constraint.isFulfilled(point)`	whether L ≤ f(x) ≤ U
`calculateGradient(point)` / `calculateHessian(point)`	the nonzeros as `{index: value}` / `{(index1, index2): value}`
`expression.calculate(point)`	the value of an expression, also one that is not in a constraint
`problem.getMostDeviatingNumericConstraint(point)`	the `NumericConstraintValue` of the most violated constraint, or `None`

A point without a value for every variable raises ValueError.

The Hessian is symmetric, so calculateHessian only returns the upper triangle, i.e., the elements with index1 ≤ index2. The example below mirrors them into the full matrix.

problem = SHOTpy.Problem(solver)

x = problem.addVariable("x", SHOTpy.VariableType.Real, 0.0, 10.0)
y = problem.addVariable("y", SHOTpy.VariableType.Real, 0.1, 10.0)

# 1 <= x + 2*y <= 8
problem.addConstraint(SHOTpy.inequality(1, x + 2 * y, 8), "lin")

# exp(x) + log(y) <= 5
problem.addConstraint(SHOTpy.exp(x) + SHOTpy.log(y) <= 5, "nl")

# minimize x^2 - y + 4
problem.setObjective(x**2 - y + 4, SHOTpy.ObjectiveDirection.Minimize)
problem.finalize()

# finalize() can replace or rewrite the constraints and the objective, so they are read back from the problem
constraints = {c.name: c for c in problem.numericConstraints}
objective = problem.objectiveFunction

point = [2.0, 1.0]
print(f"Objective value   : {objective.calculateValue(point):.4f}")
print(f"Objective gradient: {objective.calculateGradient(point)}")

for name, constraint in constraints.items():
    value = constraint.calculateNumericValue(point)
    print(f"{constraint.name}: f(x) = {value.functionValue:.4f}, fulfilled = {value.isFulfilled}, "
          f"error = {value.error:.4f}, gradient = {constraint.calculateGradient(point)}")

worst = problem.getMostDeviatingNumericConstraint(point)
print(f"Most violated constraint: {worst.constraint.name if worst else None}")

# The Hessian only has the upper triangle, so it is mirrored to get the full matrix
hessian = constraints["nl"].calculateHessian(point)
print(f"Hessian of nl     : {hessian}")

matrix = [[0.0] * len(point) for _ in point]
for (i, j), value in hessian.items():
    matrix[i][j] = matrix[j][i] = value
print(f"As a full matrix  : {matrix}")

# An expression can be evaluated directly, without adding it to a constraint
print(f"x*y + exp(x)      : {(x * y + SHOTpy.exp(x)).calculate(point):.4f}")

Interval Evaluation¶

The functions can also be bounded over a box, which is a list with an Interval for each variable. A tuple (lower, upper) or a number can be used in place of an Interval:

Method	Returns
`objectiveFunction.calculateValue(box)` / `objectiveFunction.getBounds()`	the bounds over the box / over the variable bounds
`constraint.calculateFunctionValue(box)` / `constraint.getConstraintFunctionBounds()`	the bounds of f(x) over the box / over the variable bounds
`expression.calculate(box)` / `expression.getBounds()`	the bounds of an expression over the box / over the variable bounds
`problem.getVariableBounds()` / `variable.getBound()`	the variable bounds as `Interval`s

The bounds are calculated with interval arithmetic, so they contain all the values of the function in the box, but are often wider than the range of the values. The arguments of functions such as log and sqrt are restricted to their domains. The bounds of a constraint or the objective function are [SHOT_DBL_MIN, SHOT_DBL_MAX] if they cannot be calculated, while an expression raises SHOTpy.IntervalError. The bounds are not rounded outwards, so they are not verified enclosures of the floating-point calculations.

box = [(1.0, 3.0), SHOTpy.Interval(0.5, 2.0)]

print(f"Objective bounds over the box   : {objective.calculateValue(box)}")
print(f"Objective bounds over the bounds: {objective.getBounds()}")

for name, constraint in constraints.items():
    print(f"{name}: f(x) in {constraint.calculateFunctionValue(box)} over the box")

# An Interval supports arithmetic and can be unpacked
bounds = (x * y + SHOTpy.exp(x)).calculate(box)
lower, upper = bounds
print(f"x*y + exp(x) in [{lower:.4f}, {upper:.4f}], width {bounds.width:.4f}, contains 10: {10 in bounds}")
print(f"Twice the bounds: {2 * bounds}")

try:
    SHOTpy.log(y - 5).calculate(box)
except SHOTpy.IntervalError as error:
    print(f"IntervalError: {error}")

Interval Derivatives and a Numerical Convexity Check¶

calculateIntervalDerivatives(box, order=2) of a constraint, the objective function or an expression returns an IntervalDerivatives with bounds over the box of the value (value), the gradient (gradient, {index: Interval}) and the upper triangle of the Hessian (hessian, {(index1, index2): Interval}). The bounds only enclose the derivatives when the function is defined everywhere in the box (isDefinedOnBox) and continuously differentiable up to the order of the derivative (differentiableOrder); isBounded tells whether all the bounds are finite.

A twice continuously differentiable function is convex on a box if its Hessian is positive semidefinite everywhere in it. The interval Hessian contains all the Hessians in the box, so by Gershgorin’s theorem it is enough that each diagonal element is at least the sum of the absolute values of the other elements of its row:

\[ \operatorname{lower}(H_{ii}) - \sum_{j\ne i}\max\!\left(\left|\operatorname{lower}(H_{ij})\right|,\left|\operatorname{upper}(H_{ij})\right|\right) \ge 0 \quad \text{for every }i. \]

The bounds are not rounded outwards, so this is a numerical check, not a proof: a verified certificate needs outward rounding both in SHOT’s interval arithmetic and in the sums of the check.

def check_convexity(derivatives, number_of_variables):
    """True if the function passes the convexity check on the box, False if it does not, which is inconclusive, and
    None if the bounds cannot be used, e.g., since the function is not twice continuously differentiable there."""
    if not (derivatives.isDefinedOnBox and derivatives.differentiableOrder == 2 and derivatives.isBounded):
        return None

    hessian = derivatives.hessian

    for i in range(number_of_variables):
        diagonal = hessian[(i, i)].lower if (i, i) in hessian else 0.0
        others = sum(max(abs(h.lower), abs(h.upper)) for (j, k), h in hessian.items() if (j == i) != (k == i))

        if diagonal - others < 0.0:
            return False

    return True

# x*y is not convex, so SHOT classifies exp(x) + exp(y) + x*y as nonconvex, which it is on, e.g., [-5, -4]^2.
# On [1, 2]^2, however, the diagonal of the Hessian [[exp(x), 1], [1, exp(y)]] is at least e > 1.
function = SHOTpy.exp(x) + SHOTpy.exp(y) + x * y
derivatives = function.calculateIntervalDerivatives([(1, 2), (1, 2)])
print(derivatives)
print(f"Gradient: {derivatives.gradient}")
print(f"Hessian : {derivatives.hessian}")
print(f"Convex on [1, 2]^2  : {check_convexity(derivatives, 2)}")
print(f"Convex on [-5, -4]^2: {check_convexity(function.calculateIntervalDerivatives([(-5, -4), (-5, -4)]), 2)}")

# x^3 + y^2 is convex where x >= 0. Splitting the box shows where the check passes
function = x**3 + y**2

for k in range(6):
    box = [(-1.0 + 0.5 * k, -0.5 + 0.5 * k), (-1.0, 1.0)]
    print(f"x in [{box[0][0]:4.1f}, {box[0][1]:4.1f}]: {check_convexity(function.calculateIntervalDerivatives(box), 2)}")

# |x| + y^2 is convex, but not twice differentiable at x = 0, so the bounds cannot be used there
derivatives = (SHOTpy.abs(x) + y**2).calculateIntervalDerivatives([(-1, 1), (-1, 1)])
print(f"|x| + y^2: differentiable order {derivatives.differentiableOrder}, check {check_convexity(derivatives, 2)}")

import math

# A failed check is inconclusive, since the bounds can be wider than the Hessians. A point in the interior of the box
# where the Hessian has a negative eigenvalue shows that the function is not convex on the box.
def eigenvalues(hessian):
    """The eigenvalues of a symmetric 2x2 matrix given by its upper triangle."""
    a, b, c = hessian.get((0, 0), 0.0), hessian.get((0, 1), 0.0), hessian.get((1, 1), 0.0)
    radius = math.sqrt(((a - c) / 2) ** 2 + b**2)
    return ((a + c) / 2 - radius, (a + c) / 2 + radius)

bilinear_problem = SHOTpy.Problem(solver)
u = bilinear_problem.addVariable("u", SHOTpy.VariableType.Real, -1.0, 1.0)
v = bilinear_problem.addVariable("v", SHOTpy.VariableType.Real, -1.0, 1.0)
bilinear_problem.addConstraint(u * v <= 0.5, "bilinear")
bilinear_problem.setObjective(u + v, SHOTpy.ObjectiveDirection.Minimize)
bilinear_problem.finalize()

bilinear = bilinear_problem.numericConstraints[0]
print(f"Check on [-1, 1]^2         : {check_convexity(bilinear.calculateIntervalDerivatives([(-1, 1), (-1, 1)]), 2)}")
print(f"Eigenvalues at (0.2, -0.3): {eigenvalues(bilinear.calculateHessian([0.2, -0.3]))}")

αBB Underestimators¶

The interval Hessian also gives a convex underestimator of a nonconvex function on a box, as in the αBB method. For a twice continuously differentiable f on the box [l, u],

\[ L(x) = f(x) + \alpha \sum_i (x_i-l_i)(x_i-u_i) \]

is below f on the box, since each added term is at most zero, and equal to f at the corners of the box. Its Hessian is that of f plus 2αI, so L is convex if 2α is at least minus the smallest eigenvalue of the Hessian anywhere in the box. The Gershgorin bound of the convexity check gives a lower bound on that eigenvalue:

\[ \alpha = \max\!\left(0,-\frac12\min_i\left[\operatorname{lower}(H_{ii})-\sum_{j\ne i}\max\!\left(\left|\operatorname{lower}(H_{ij})\right|,\left|\operatorname{upper}(H_{ij})\right|\right)\right]\right). \]

For a function of one variable, this is \(\alpha=\max(0,-\tfrac12\operatorname{lower}(f^{\prime\prime}))\). The largest gap between f and L is \(\alpha(u-l)^2/4\), in the middle of the interval, so it shrinks quadratically when the box is split, which is what a branch-and-bound method using αBB relies on. As for the convexity check, the bounds are not rounded outwards, so L is only numerically an underestimator.

def alpha_bb(derivatives, number_of_variables):
    """The alpha of the alphaBB underestimator from the Gershgorin bound of the interval Hessian."""
    assert derivatives.isDefinedOnBox and derivatives.differentiableOrder == 2 and derivatives.isBounded
    hessian = derivatives.hessian
    smallest = min((hessian[(i, i)].lower if (i, i) in hessian else 0.0)
                   - sum(max(abs(h.lower), abs(h.upper)) for (j, k), h in hessian.items() if (j == i) != (k == i))
                   for i in range(number_of_variables))
    return max(0.0, -0.5 * smallest)

abb_problem = SHOTpy.Problem(solver)
t = abb_problem.addVariable("t", SHOTpy.VariableType.Real, -2.0, 2.0)
f = SHOTpy.sin(3 * t) + 0.25 * t**2   # nonconvex: f''(t) = -9 sin(3t) + 0.5

def underestimator(lower, upper):
    """alpha and L(t) = f(t) + alpha (t - lower)(t - upper) on [lower, upper]."""
    alpha = alpha_bb(f.calculateIntervalDerivatives([(lower, upper)]), 1)
    return alpha, (lambda value: f.calculate([value]) + alpha * (value - lower) * (value - upper))

pieces = {"[-2, 2]": [(-2.0, 2.0)], "[-2, 0] and [0, 2]": [(-2.0, 0.0), (0.0, 2.0)]}
points = [-2.0 + 4.0 * k / 400 for k in range(401)]

for name, boxes in pieces.items():
    largest_gap = 0.0
    for lower, upper in boxes:
        alpha, L = underestimator(lower, upper)
        inside = [p for p in points if lower <= p <= upper]
        gaps = [f.calculate([p]) - L(p) for p in inside]
        assert min(gaps) >= -1e-12                                    # L is below f
        assert abs(L(lower) - f.calculate([lower])) < 1e-12          # and equal to f at the ends
        largest_gap = max(largest_gap, max(gaps))
        print(f"[{lower:4.1f}, {upper:4.1f}]: alpha = {alpha:.4f}, largest gap {max(gaps):.4f}"
              f" (alpha (u - l)^2 / 4 = {alpha * (upper - lower) ** 2 / 4:.4f})")

    print(f"Underestimator on {name}: largest gap {largest_gap:.4f}")

# The underestimator on [-2, 2] passes the numerical convexity check
alpha, _ = underestimator(-2.0, 2.0)
convexified = f + alpha * (t + 2.0) * (t - 2.0)
assert check_convexity(convexified.calculateIntervalDerivatives([(-2.0, 2.0)]), 1)

try:
    import matplotlib.pyplot as plt
except ImportError:
    plt = None

if plt is not None:
    figure, axes = plt.subplots(figsize=(8, 4.5))
    axes.plot(points, [f.calculate([p]) for p in points], color="black", linewidth=2, label="f(t) = sin(3t) + t²/4")

    for (name, boxes), color in zip(pieces.items(), ["tab:blue", "tab:orange"]):
        for number, (lower, upper) in enumerate(boxes):
            _, L = underestimator(lower, upper)
            inside = [p for p in points if lower <= p <= upper]
            axes.plot(inside, [L(p) for p in inside], color=color, linestyle="--",
                      label=f"αBB underestimator on {name}" if number == 0 else None)

    axes.axvline(0.0, color="gray", linewidth=0.5)
    axes.set_xlabel("t")
    axes.set_title("αBB underestimators from the interval Hessian")
    axes.legend()
    plt.show()