Unexpected expression type error when using new piecewise linear function

[schmidtmj]

Summary

Hi. I'm new to Pyomo so this may be a user error or a true bug.

I was working with the new pyomo.contrib piecewise library with the hope of using its multi-dimensional piecewise capabilities for an upcoming project. I am able to generate a piecewise function and then use this function to generate an indexed constraint but get an error when trying to solve the model.

It appears that the PiecewiseLinearExpression is an unexpected expression when pyomo generates a standard representation before executing the solver. I got the same error using the 'GLPK' or 'GUROBI' solvers. Maybe this is an issue in how the piecewise linear expression is being generated or registered? but I'm not really a developer so can't quite sort out what has gone wrong.

Happy to answer any other questions. Thanks.

See minimal reproducible example below.

Steps to reproduce the issue

import pyomo.environ, pyomo.core as pyo
from pyomo.contrib.piecewise import PiecewiseLinearFunction

def f(x):
    y = x**2
    return y

def r_1(model, t):
    return model.Y[t] == model.f_1(model.X[t])

def r_obj(model):
    return sum([model.Y[t] for t in model.sT])

def create_model():

    t_list = list(range(10))
    x_max = 100
    x_min = 0
    pw_pts = [x_min, 10, 50, x_max]

    model = pyo.ConcreteModel()
    model.sT = pyo.Set(initialize=t_list)
    model.X = pyo.Var(model.sT, domain=pyo.NonNegativeReals, bounds=(0, x_max))
    model.Y = pyo.Var(model.sT, domain=pyo.NonNegativeReals)

    model.f_1 = PiecewiseLinearFunction(points=pw_pts, function=f)
    model.c_1 = pyo.Constraint(model.sT, expr=r_1)
    model.obj = pyo.Objective(rule=r_obj, sense=pyo.maximize)

    return model

model = create_model()
solver = pyomo.environ.SolverFactory('gurobi')
solution = solver.solve(model)

Error Message

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
Cell In[5], line 35
     33 model = create_model()
     34 solver = pyomo.environ.SolverFactory('gurobi')
---> 35 solution = solver.solve(model)

File ~/miniconda3/envs/pyomo/lib/python3.10/site-packages/pyomo/opt/base/solvers.py:597, in OptSolver.solve(self, *args, **kwds)
    593 try:
    594     # we're good to go.
    595     initial_time = time.time()
--> 597     self._presolve(*args, **kwds)
    599     presolve_completion_time = time.time()
    600     if self._report_timing:

File ~/miniconda3/envs/pyomo/lib/python3.10/site-packages/pyomo/solvers/plugins/solvers/GUROBI.py:230, in GUROBISHELL._presolve(self, *args, **kwds)
    225         self._warm_start_file_name = TempfileManager.create_tempfile(
    226             suffix='.gurobi.mst'
    227         )
    229 # let the base class handle any remaining keywords/actions.
--> 230 ILMLicensedSystemCallSolver._presolve(self, *args, **kwds)
    232 # NB: we must let the base class presolve run first so that the
    233 # symbol_map is actually constructed!
    235 if (len(args) > 0) and (not isinstance(args[0], str)):

File ~/miniconda3/envs/pyomo/lib/python3.10/site-packages/pyomo/opt/solver/shellcmd.py:224, in SystemCallSolver._presolve(self, *args, **kwds)
    221 self._keepfiles = kwds.pop("keepfiles", False)
    222 self._define_signal_handlers = kwds.pop('use_signal_handling', None)
--> 224 OptSolver._presolve(self, *args, **kwds)
    226 #
    227 # Verify that the input problems exists
    228 #
    229 for filename in self._problem_files:

File ~/miniconda3/envs/pyomo/lib/python3.10/site-packages/pyomo/opt/base/solvers.py:706, in OptSolver._presolve(self, *args, **kwds)
    700 if self._problem_format:
    701     write_start_time = time.time()
    702     (
    703         self._problem_files,
    704         self._problem_format,
    705         self._smap_id,
--> 706     ) = self._convert_problem(
    707         args, self._problem_format, self._valid_problem_formats, **kwds
    708     )
    709     total_time = time.time() - write_start_time
    710     if self._report_timing:

File ~/miniconda3/envs/pyomo/lib/python3.10/site-packages/pyomo/opt/base/solvers.py:757, in OptSolver._convert_problem(self, args, problem_format, valid_problem_formats, **kwds)
    756 def _convert_problem(self, args, problem_format, valid_problem_formats, **kwds):
--> 757     return convert_problem(
    758         args, problem_format, valid_problem_formats, self.has_capability, **kwds
    759     )

File ~/miniconda3/envs/pyomo/lib/python3.10/site-packages/pyomo/opt/base/convert.py:99, in convert_problem(args, target_problem_type, valid_problem_types, has_capability, **kwds)
     97                 tmpkw = kwds
     98                 tmpkw['capabilities'] = has_capability
---> 99                 problem_files, symbol_map = converter.apply(*tmp, **tmpkw)
    100                 return problem_files, ptype, symbol_map
    102 msg = 'No conversion possible.  Source problem type: %s.  Valid target types: %s'

File ~/miniconda3/envs/pyomo/lib/python3.10/site-packages/pyomo/solvers/plugins/converter/model.py:78, in PyomoMIPConverter.apply(self, *args, **kwds)
     70         symbol_map_id = instance.write(
     71             problem_filename,
     72             format=ProblemFormat.cpxlp,
   (...)
     75             **io_options
     76         )
     77     else:
---> 78         (problem_filename, symbol_map_id) = instance.write(
     79             filename=problem_filename,
     80             format=ProblemFormat.cpxlp,
     81             solver_capability=capabilities,
     82             io_options=io_options,
     83         )
     84     return (problem_filename,), symbol_map_id
     85 else:
     86     #
     87     # I'm simply exposing a fatal issue with
   (...)
     90     # arguments that can be sent to the writer?
     91     #

File ~/miniconda3/envs/pyomo/lib/python3.10/site-packages/pyomo/core/base/block.py:2086, in _BlockData.write(self, filename, format, solver_capability, io_options)
   2083     def solver_capability(x):
   2084         return True
-> 2086 (filename, smap) = problem_writer(self, filename, solver_capability, io_options)
   2087 smap_id = id(smap)
   2088 if not hasattr(self, 'solutions'):
   2089     # This is a bit of a hack.  The write() method was moved
   2090     # here from PyomoModel to support the solution of arbitrary
   (...)
   2095     # dependency (we only need it here because we store the
   2096     # SymbolMap returned by the writer in the solutions).

File ~/miniconda3/envs/pyomo/lib/python3.10/site-packages/pyomo/repn/plugins/cpxlp.py:164, in ProblemWriter_cpxlp.__call__(self, model, output_filename, solver_capability, io_options)
    162 with PauseGC() as pgc:
    163     with open(output_filename, "w") as output_file:
--> 164         symbol_map = self._print_model_LP(
    165             model,
    166             output_file,
    167             solver_capability,
    168             labeler,
    169             output_fixed_variable_bounds=output_fixed_variable_bounds,
    170             file_determinism=file_determinism,
    171             row_order=row_order,
    172             column_order=column_order,
    173             skip_trivial_constraints=skip_trivial_constraints,
    174             force_objective_constant=force_objective_constant,
    175             include_all_variable_bounds=include_all_variable_bounds,
    176         )
    178 self._referenced_variable_ids.clear()
    180 return output_filename, symbol_map

File ~/miniconda3/envs/pyomo/lib/python3.10/site-packages/pyomo/repn/plugins/cpxlp.py:700, in ProblemWriter_cpxlp._print_model_LP(self, model, output_file, solver_capability, labeler, output_fixed_variable_bounds, file_determinism, row_order, column_order, skip_trivial_constraints, force_objective_constant, include_all_variable_bounds)
    697     yield_all_constraints = constraint_generator
    699 # FIXME: This is a hack to get nested blocks working...
--> 700 for constraint_data, repn in yield_all_constraints():
    701     have_nontrivial = True
    703     degree = repn.polynomial_degree()

File ~/miniconda3/envs/pyomo/lib/python3.10/site-packages/pyomo/repn/plugins/cpxlp.py:680, in ProblemWriter_cpxlp._print_model_LP.<locals>.constraint_generator()
    678     repn = constraint_data.canonical_form()
    679 else:
--> 680     repn = generate_standard_repn(
    681         constraint_data.body,
    682         quadratic=supports_quadratic_constraint,
    683     )
    684 if gen_con_repn:
    685     block_repn[constraint_data] = repn

File ~/miniconda3/envs/pyomo/lib/python3.10/site-packages/pyomo/repn/standard_repn.py:384, in generate_standard_repn(expr, idMap, compute_values, verbose, quadratic, repn)
    368     raise ValueError("Unexpected expression type: " + str(expr))
    370 #
    371 # WEH - Checking the polynomial degree didn't
    372 #       turn out to be a win.  But I'm leaving this
   (...)
    382 #                        repn=repn)
    383 # else:
--> 384 return _generate_standard_repn(
    385     expr,
    386     idMap=idMap,
    387     compute_values=compute_values,
    388     verbose=verbose,
    389     quadratic=quadratic,
    390     repn=repn,
    391 )

File ~/miniconda3/envs/pyomo/lib/python3.10/site-packages/pyomo/repn/standard_repn.py:1180, in _generate_standard_repn(expr, idMap, compute_values, verbose, quadratic, repn)
   1173 def _generate_standard_repn(
   1174     expr, idMap=None, compute_values=True, verbose=False, quadratic=True, repn=None
   1175 ):
   1176     if expr.__class__ is EXPR.SumExpression:
   1177         #
   1178         # This is the common case, so start collecting the sum
   1179         #
-> 1180         ans = _collect_sum(expr, 1, idMap, compute_values, verbose, quadratic)
   1181     else:
   1182         #
   1183         # Call generic recursive logic
   1184         #
   1185         ans = _collect_standard_repn(expr, 1, idMap, compute_values, verbose, quadratic)

File ~/miniconda3/envs/pyomo/lib/python3.10/site-packages/pyomo/repn/standard_repn.py:509, in _collect_sum(exp, multiplier, idMap, compute_values, verbose, quadratic)
    507         ans.constant += multiplier * e_
    508 else:
--> 509     res_ = _collect_standard_repn(
    510         e_, multiplier, idMap, compute_values, verbose, quadratic
    511     )
    512     #
    513     # Add returned from recursion
    514     #
    515     ans.constant += res_.constant

File ~/miniconda3/envs/pyomo/lib/python3.10/site-packages/pyomo/repn/standard_repn.py:1147, in _collect_standard_repn(exp, multiplier, idMap, compute_values, verbose, quadratic)
   1145 fn = _repn_collectors.get(exp.__class__, None)
   1146 if fn is not None:
-> 1147     return fn(exp, multiplier, idMap, compute_values, verbose, quadratic)
   1148 #
   1149 # Catch any known numeric constants
   1150 #
   1151 if exp.__class__ in native_numeric_types or not exp.is_potentially_variable():

File ~/miniconda3/envs/pyomo/lib/python3.10/site-packages/pyomo/repn/standard_repn.py:1022, in _collect_negation(exp, multiplier, idMap, compute_values, verbose, quadratic)
   1021 def _collect_negation(exp, multiplier, idMap, compute_values, verbose, quadratic):
-> 1022     return _collect_standard_repn(
   1023         exp._args_[0], -1 * multiplier, idMap, compute_values, verbose, quadratic
   1024     )

File ~/miniconda3/envs/pyomo/lib/python3.10/site-packages/pyomo/repn/standard_repn.py:1147, in _collect_standard_repn(exp, multiplier, idMap, compute_values, verbose, quadratic)
   1145 fn = _repn_collectors.get(exp.__class__, None)
   1146 if fn is not None:
-> 1147     return fn(exp, multiplier, idMap, compute_values, verbose, quadratic)
   1148 #
   1149 # Catch any known numeric constants
   1150 #
   1151 if exp.__class__ in native_numeric_types or not exp.is_potentially_variable():

File ~/miniconda3/envs/pyomo/lib/python3.10/site-packages/pyomo/repn/standard_repn.py:1045, in _collect_identity(exp, multiplier, idMap, compute_values, verbose, quadratic)
   1043     else:
   1044         return Results(constant=multiplier * exp._args_[0])
-> 1045 return _collect_standard_repn(
   1046     exp.expr, multiplier, idMap, compute_values, verbose, quadratic
   1047 )

File ~/miniconda3/envs/pyomo/lib/python3.10/site-packages/pyomo/repn/standard_repn.py:1168, in _collect_standard_repn(exp, multiplier, idMap, compute_values, verbose, quadratic)
   1166     _repn_collectors[exp.__class__] = fn
   1167     return fn(exp, multiplier, idMap, compute_values, verbose, quadratic)
-> 1168 raise ValueError(
   1169     "Unexpected expression (type %s)" % type(exp).__name__
   1170 )

ValueError: Unexpected expression (type PiecewiseLinearExpression)

Information on your system

Pyomo version: 6.5.1a0 Python version: 3.10.11, running in Jupyter Notebook Operating system: Ubuntu 18.04.6 How Pyomo was installed (PyPI, conda, source): Conda environment with Pip to get latest Pyomo Solver (if applicable): Tried with Gurobi and GLPK and got same result.

Additional information

[jsiirola]

This is not a bug (although the error you get from generate_standard_repn() is admittedly less than useful -- that will be at least somewhat resolved by #2823).

The short answer is that you have to transform your model into a MIP before sending it to a MIP solver. There are multiple ways (formulations) that you could use when represending the problem as a MIP, and Pyomo won't guess which one is most appropriate for your model. The new piecewise library provides utilities to convert the model to a MIP in two steps: first a transformation to replace the PiecewiseLinearFunction with equivalent Disjunctions (now the model is a GDP), then in a second step, the GDP is converted to a MIP.

This package is still under development, and we haven't had a chance to complete the documentation for it yet. Until then, you can see the set of supported transformations in pyomo/contrib/piecewise/transform. This includes several conversions from piecewise to GDP (including contrib.piecewise.outer_repn_gdp, contrib.piecewise.inner_repn_gdp, and contrib.piecewise.reduced_inner_repn_gdp), as well as several "convenience" transformations that call both a Piecewise-to-GDP transformation and a GDP-to-MIP transformation (e.g., contrib.convex_combination, contrib.piecewise.disaggregated_convex_combination, and contrib.piecewise.multiple_choice)

[schmidtmj]

Thanks for the kind, informative, and speedy response.