Expressions

In this section, we use the word “expression” in two ways: first in the general sense of the word and second to describe a class of Pyomo objects that have the name Expression as described in the subsection on expression objects.

Rules to Generate Expressions

Both objectives and constraints make use of rules to generate expressions. These are Python functions that return the appropriate expression. These are first-class functions that can access global data as well as data passed in, including the model object.

Operations on model elements results in expressions, which seems natural in expressions like the constraints we have seen so far. It is also possible to build up expressions. The following example illustrates this, along with a reference to global Python data in the form of a Python variable called switch:

switch = 3

model.A = RangeSet(1, 10)
model.c = Param(model.A)
model.d = Param()
model.x = Var(model.A, domain=Boolean)


def pi_rule(model):
    accexpr = summation(model.c, model.x)
    if switch >= 2:
        accexpr = accexpr - model.d
    return accexpr >= 0.5


PieSlice = Constraint(rule=pi_rule)

In this example, the constraint that is generated depends on the value of the Python variable called switch. If the value is 2 or greater, then the constraint is summation(model.c, model.x) - model.d >= 0.5; otherwise, the model.d term is not present.

Warning

Because model elements result in expressions, not values, the following does not work as expected in an abstract model!

model.A = RangeSet(1, 10)
model.c = Param(model.A)
model.d = Param()
model.x = Var(model.A, domain=Boolean)


def pi_rule(model):
    accexpr = summation(model.c, model.x)
    if model.d >= 2:  # NOT in an abstract model!!
        accexpr = accexpr - model.d
    return accexpr >= 0.5


PieSlice = Constraint(rule=pi_rule)

The trouble is that model.d >= 2 results in an expression, not its evaluated value. Instead use if value(model.d) >= 2

Note

Pyomo supports non-linear expressions and can call non-linear solvers such as Ipopt.

Piecewise Linear Expressions

Pyomo has facilities to add piecewise constraints of the form y=f(x) for a variety of forms of the function f.

The piecewise types other than SOS2, BIGM_SOS1, BIGM_BIN are implement as described in the paper [Vielma_et_al].

There are two basic forms for the declaration of the constraint:

# model.pwconst = Piecewise(indexes, yvar, xvar, **Keywords)
# model.pwconst = Piecewise(yvar,xvar,**Keywords)

where pwconst can be replaced by a name appropriate for the application. The choice depends on whether the x and y variables are indexed. If so, they must have the same index sets and these sets are give as the first arguments.

Keywords:

  • pw_pts={ },[ ],( )

    A dictionary of lists (where keys are the index set) or a single list (for the non-indexed case or when an identical set of breakpoints is used across all indices) defining the set of domain breakpoints for the piecewise linear function.

    Note

    pw_pts is always required. These give the breakpoints for the piecewise function and are expected to fully span the bounds for the independent variable(s).

  • pw_repn=<Option>

    Indicates the type of piecewise representation to use. This can have a major impact on solver performance. Options: (Default “SOS2”)

    • “SOS2” - Standard representation using sos2 constraints.

    • “BIGM_BIN” - BigM constraints with binary variables. The theoretically tightest M values are automatically determined.

    • “BIGM_SOS1” - BigM constraints with sos1 variables. The theoretically tightest M values are automatically determined.

    • “DCC” - Disaggregated convex combination model.

    • “DLOG” - Logarithmic disaggregated convex combination model.

    • “CC” - Convex combination model.

    • “LOG” - Logarithmic branching convex combination.

    • “MC” - Multiple choice model.

    • “INC” - Incremental (delta) method.

    Note

    Step functions are supported for all but the two BIGM options. Refer to the ‘force_pw’ option.

  • pw_constr_type= <Option>

    Indicates the bound type of the piecewise function. Options:

    • “UB” - y variable is bounded above by piecewise function.

    • “LB” - y variable is bounded below by piecewise function.

    • “EQ” - y variable is equal to the piecewise function.

  • f_rule=f(model,i,j,…,x), { }, [ ], ( )

    An object that returns a numeric value that is the range value corresponding to each piecewise domain point. For functions, the first argument must be a Pyomo model. The last argument is the domain value at which the function evaluates (Not a Pyomo Var). Intermediate arguments are the corresponding indices of the Piecewise component (if any). Otherwise, the object can be a dictionary of lists/tuples (with keys the same as the indexing set) or a singe list/tuple (when no indexing set is used or when all indices use an identical piecewise function). Examples:

    # A function that changes with index
    def f(model, j, x):
        if j == 2:
            return x**2 + 1.0
        else:
            return x**2 + 5.0
    
    
    # A nonlinear function
    f = lambda model, x: exp(x) + value(model.p)
    
    # A step function
    f = [0, 0, 1, 1, 2, 2]
    
  • force_pw=True/False

    Using the given function rule and pw_pts, a check for convexity/concavity is implemented. If (1) the function is convex and the piecewise constraints are lower bounds or if (2) the function is concave and the piecewise constraints are upper bounds then the piecewise constraints will be substituted for linear constraints. Setting ‘force_pw=True’ will force the use of the original piecewise constraints even when one of these two cases applies.

  • warning_tol=<float>

    To aid in debugging, a warning is printed when consecutive slopes of piecewise segments are within <warning_tol> of each other. Default=1e-8

  • warn_domain_coverage=True/False

    Print a warning when the feasible region of the domain variable is not completely covered by the piecewise breakpoints. Default=True

  • unbounded_domain_var=True/False

    Allow an unbounded or partially bounded Pyomo Var to be used as the domain variable. Default=False

    Note

    This does not imply unbounded piecewise segments will be constructed. The outermost piecewise breakpoints will bound the domain variable at each index. However, the Var attributes .lb and .ub will not be modified.

Here is an example of an assignment to a Python dictionary variable that has keywords for a picewise constraint:

kwds = {'pw_constr_type': 'EQ', 'pw_repn': 'SOS2', 'sense': maximize, 'force_pw': True}

Here is a simple example based on the example given earlier in Symbolic Index Sets. In this new example, the objective function is the sum of c times x to the fourth. In this example, the keywords are passed directly to the Piecewise function without being assigned to a dictionary variable. The upper bound on the x variables was chosen whimsically just to make the example. The important thing to note is that variables that are going to appear as the independent variable in a piecewise constraint must have bounds.

#  ___________________________________________________________________________
#
#  Pyomo: Python Optimization Modeling Objects
#  Copyright (c) 2008-2024
#  National Technology and Engineering Solutions of Sandia, LLC
#  Under the terms of Contract DE-NA0003525 with National Technology and
#  Engineering Solutions of Sandia, LLC, the U.S. Government retains certain
#  rights in this software.
#  This software is distributed under the 3-clause BSD License.
#  ___________________________________________________________________________

# abstract2piece.py
# Similar to abstract2.py, but the objective is now c times x to the fourth power

from pyomo.environ import *

model = AbstractModel()

model.I = Set()
model.J = Set()

Topx = 6.1  # range of x variables

model.a = Param(model.I, model.J)
model.b = Param(model.I)
model.c = Param(model.J)

# the next line declares a variable indexed by the set J
model.x = Var(model.J, domain=NonNegativeReals, bounds=(0, Topx))
model.y = Var(model.J, domain=NonNegativeReals)

# to avoid warnings, we set breakpoints at or beyond the bounds
PieceCnt = 100
bpts = []
for i in range(PieceCnt + 2):
    bpts.append(float((i * Topx) / PieceCnt))


def f4(model, j, xp):
    # we not need j, but it is passed as the index for the constraint
    return xp**4


model.ComputeObj = Piecewise(
    model.J, model.y, model.x, pw_pts=bpts, pw_constr_type='EQ', f_rule=f4
)


def obj_expression(model):
    return summation(model.c, model.y)


model.OBJ = Objective(rule=obj_expression)


def ax_constraint_rule(model, i):
    # return the expression for the constraint for i
    return sum(model.a[i, j] * model.x[j] for j in model.J) >= model.b[i]


# the next line creates one constraint for each member of the set model.I
model.AxbConstraint = Constraint(model.I, rule=ax_constraint_rule)

A more advanced example is provided in abstract2piecebuild.py in BuildAction and BuildCheck.

Expression Objects

Pyomo Expression objects are very similar to the Param component (with mutable=True) except that the underlying values can be numeric constants or Pyomo expressions. Here’s an illustration of expression objects in an AbstractModel. An expression object with an index set that is the numbers 1, 2, 3 is created and initialized to be the model variable x times the index. Later in the model file, just to illustrate how to do it, the expression is changed but just for the first index to be x squared.

model = ConcreteModel()
model.x = Var(initialize=1.0)


def _e(m, i):
    return m.x * i


model.e = Expression([1, 2, 3], rule=_e)

instance = model.create_instance()

print(value(instance.e[1]))  # -> 1.0
print(instance.e[1]())  # -> 1.0
print(instance.e[1].value)  # -> a pyomo expression object

# Change the underlying expression
instance.e[1].value = instance.x**2

# ... solve
# ... load results

# print the value of the expression given the loaded optimal solution
print(value(instance.e[1]))

An alternative is to create Python functions that, potentially, manipulate model objects. E.g., if you define a function

def f(x, p):
    return x + p


You can call this function with or without Pyomo modeling components as the arguments. E.g., f(2,3) will return a number, whereas f(model.x, 3) will return a Pyomo expression due to operator overloading.

If you take this approach you should note that anywhere a Pyomo expression is used to generate another expression (e.g., f(model.x, 3) + 5), the initial expression is always cloned so that the new generated expression is independent of the old. For example:

model = ConcreteModel()
model.x = Var()

# create a Pyomo expression
e1 = model.x + 5

# create another Pyomo expression
# e1 is copied when generating e2
e2 = e1 + model.x

If you want to create an expression that is shared between other expressions, you can use the Expression component.