nonlinear_to_pwl triangulation fails on some models #3763
Description
Summary
For some minlplib problems, the nonlinear_to_pwl transform fails during the Delaunay, J1, and OrderedJ1 triangulation methods.
Steps to reproduce the issue
As an example, take the maxmin problem. (I attached a directory at the bottom with the problem and scripts to reproduce).
from pathlib import Path
import pickle
import pyomo.environ as pe
from pyomo.contrib import piecewise as pw
def load_model(case_name):
pickle_path = Path.cwd().absolute().joinpath(f"{case_name}.pkl")
print(pickle_path)
assert pickle_path.exists()
assert pickle_path.is_file()
with open(pickle_path, 'rb') as f:
m = pickle.load(f)
return m
def pw_trans(m, num_points=3):
trans = pe.TransformationFactory("contrib.piecewise.nonlinear_to_pwl")
trans.apply_to(
m,
num_points=num_points,
targets=list(m.component_data_objects(pe.Constraint, descend_into=True)),
additively_decompose=True,
)
if __name__ == "__main__":
m = load_model('maxmin')
pw_trans(m, num_points=3)
The triangulation is performed in piecewise_linear_function:_construct_simplices_from_multivariate_points. We can't pass the argument from nonlinear_to_pwl, but we can select different methods by adding tri = Triangulation.Delaunay # Triangulation.Delaunay | Triangulation.J1 | Triangulation.OrderedJ1 on line 333.
Error Message
At this point, we get the error:
For Triangulation.Delaunay:
ERROR: Unable to triangulate the set of input points.
ERROR: Constructing component
'_pyomo_contrib_nonlinear_to_pwl.'_pwle_cons[e1]_1'' from data=None failed:
QhullError: QH6013 qhull input error: input is less than 5-dimensional
since all points have the same x coordinate 0
While executing: | qhull d Q12 Qc Qt Qbb Qz Options selected for Qhull
2020.2.r 2020/08/31:
run-id 101538417 delaunay Q12-allow-wide Qcoplanar-keep Qtriangulate
Qbbound-last Qz-infinity-point _pre-merge _zero-centrum Qinterior-
keep Pgood _max-width 1 Error-roundoff 2.7e-15 _one-merge 3e-14
Visible-distance 1.6e-14 U-max-coplanar 1.6e-14 Width-outside 3.3e-14
_wide-facet 9.8e-14 _maxoutside 3.3e-14
Traceback (most recent call last):
File "C:/path/to/project/qhull_problem/demo_triangulation_problem_simple.py", line 31, in <module>
pw_trans(m, num_points=3)
File "C:/path/to/project/qhull_problem/demo_triangulation_problem_simple.py", line 21, in pw_trans
trans.apply_to(
File "C:/path/to/project/pyomo/pyomo/core/base/transformation.py", line 77, in apply_to
reverse_token = self._apply_to(model, **kwds)
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
File "C:/path/to/project/pyomo/pyomo/contrib/piecewise/transform/nonlinear_to_pwl.py", line 506, in _apply_to
self._apply_to_impl(instance, **kwds)
File "C:/path/to/project/pyomo/pyomo/contrib/piecewise/transform/nonlinear_to_pwl.py", line 523, in _apply_to_impl
self._transform_constraint(target, config)
File "C:/path/to/project/pyomo/pyomo/contrib/piecewise/transform/nonlinear_to_pwl.py", line 571, in _transform_constraint
pw_approx, expr_type = self._approximate_expression(
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
File "C:/path/to/project/pyomo/pyomo/contrib/piecewise/transform/nonlinear_to_pwl.py", line 691, in _approximate_expression
trans_block.add_component(f"_pwle_{name}_{k}", pwlf)
File "C:/path/to/project/pyomo/pyomo/core/base/block.py", line 1103, in add_component
val.construct(data)
File "C:/path/to/project/pyomo/pyomo/core/base/block.py", line 2245, in construct
self._getitem_when_not_present(_idx)
File "C:/path/to/project/pyomo/pyomo/contrib/piecewise/piecewise_linear_function.py", line 610, in _getitem_when_not_present
obj = handler(self, obj, parent, nonlinear_function)
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
File "C:/path/to/project/pyomo/pyomo/contrib/piecewise/piecewise_linear_function.py", line 437, in _construct_from_function_and_points
self._construct_simplices_from_multivariate_points(
File "C:/path/to/project/pyomo/pyomo/contrib/piecewise/piecewise_linear_function.py", line 337, in _construct_simplices_from_multivariate_points
triangulation = spatial.Delaunay(points)
^^^^^^^^^^^^^^^^^^^^^^^^
File "scipy/spatial/_qhull.pyx", line 1889, in scipy.spatial._qhull.Delaunay.__init__
File "scipy/spatial/_qhull.pyx", line 356, in scipy.spatial._qhull._Qhull.__init__
scipy.spatial._qhull.QhullError: QH6013 qhull input error: input is less than 5-dimensional since all points have the same x coordinate 0
While executing: | qhull d Q12 Qc Qt Qbb Qz
Options selected for Qhull 2020.2.r 2020/08/31:
run-id 101538417 delaunay Q12-allow-wide Qcoplanar-keep Qtriangulate
Qbbound-last Qz-infinity-point _pre-merge _zero-centrum Qinterior-keep
Pgood _max-width 1 Error-roundoff 2.7e-15 _one-merge 3e-14
Visible-distance 1.6e-14 U-max-coplanar 1.6e-14 Width-outside 3.3e-14
_wide-facet 9.8e-14 _maxoutside 3.3e-14
For Triangulation.J1 or Triangulation.OrderedJ1:
WARNING: LinAlgError: Singular matrix
ERROR: Constructing component
'_pyomo_contrib_nonlinear_to_pwl.'_pwle_cons[e1]_1'' from data=None failed:
DeveloperError: Internal Pyomo implementation error:
'When calculating the hyperplane approximation over the simplex
with index 0, the matrix was unexpectedly singular. This likely
means that this simplex is degenerate and that it should have been
filtered out of the triangulation'
Please report this to the Pyomo Developers.
Traceback (most recent call last):
File "C:/path/to/project/pyomo/pyomo/contrib/piecewise/piecewise_linear_function.py", line 512, in _construct_from_function_and_simplices
normal = np.linalg.solve(A, b)
^^^^^^^^^^^^^^^^^^^^^
File "C:/Python/envs/pyomo-env/Lib/site-packages/numpy/linalg/_linalg.py", line 471, in solve
r = gufunc(a, b, signature=signature)
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
File "C:/Python/envs/pyomo-env/Lib/site-packages/numpy/linalg/_linalg.py", line 163, in _raise_linalgerror_singular
raise LinAlgError("Singular matrix")
numpy.linalg.LinAlgError: Singular matrix
During handling of the above exception, another exception occurred:
Traceback (most recent call last):
File "C:/path/to/project/qhull_problem/demo_triangulation_problem_simple.py", line 30, in <module>
pw_trans(m, num_points=3)
File "C:/path/to/project/qhull_problem/demo_triangulation_problem_simple.py", line 20, in pw_trans
trans.apply_to(
File "C:/path/to/project/pyomo/pyomo/core/base/transformation.py", line 77, in apply_to
reverse_token = self._apply_to(model, **kwds)
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
File "C:/path/to/project/pyomo/pyomo/contrib/piecewise/transform/nonlinear_to_pwl.py", line 506, in _apply_to
self._apply_to_impl(instance, **kwds)
File "C:/path/to/project/pyomo/pyomo/contrib/piecewise/transform/nonlinear_to_pwl.py", line 523, in _apply_to_impl
self._transform_constraint(target, config)
File "C:/path/to/project/pyomo/pyomo/contrib/piecewise/transform/nonlinear_to_pwl.py", line 571, in _transform_constraint
pw_approx, expr_type = self._approximate_expression(
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
File "C:/path/to/project/pyomo/pyomo/contrib/piecewise/transform/nonlinear_to_pwl.py", line 691, in _approximate_expression
trans_block.add_component(f"_pwle_{name}_{k}", pwlf)
File "C:/path/to/project/pyomo/pyomo/core/base/block.py", line 1103, in add_component
val.construct(data)
File "C:/path/to/project/pyomo/pyomo/core/base/block.py", line 2245, in construct
self._getitem_when_not_present(_idx)
File "C:/path/to/project/pyomo/pyomo/contrib/piecewise/piecewise_linear_function.py", line 610, in _getitem_when_not_present
obj = handler(self, obj, parent, nonlinear_function)
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
File "C:/path/to/project/pyomo/pyomo/contrib/piecewise/piecewise_linear_function.py", line 440, in _construct_from_function_and_points
return self._construct_from_function_and_simplices(
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
File "C:/path/to/project/pyomo/pyomo/contrib/piecewise/piecewise_linear_function.py", line 525, in _construct_from_function_and_simplices
raise DeveloperError(
pyomo.common.errors.DeveloperError: Internal Pyomo implementation error:
'When calculating the hyperplane approximation over the simplex with
index 0, the matrix was unexpectedly singular. This likely means
that this simplex is degenerate and that it should have been
filtered out of the triangulation'
Please report this to the Pyomo Developers.
Information on your system
Pyomo version: 6.9.3
Python version: 3.11.13
Operating system: windows 11
How Pyomo was installed (PyPI, conda, source): source
Additional information
We can get a different error by changing the model. The following script does an expansion of the constraints before the pwl transformation:
from pathlib import Path
import pickle
import pyomo.environ as pe
from pyomo.contrib import piecewise as pw
from pyomo.core.expr.sympy_tools import sympy2pyomo_expression, sympyify_expression
from pyomo.core.expr.numeric_expr import (
NumericExpression,
SumExpression,
LinearExpression,
)
from pyomo.core.expr.numvalue import value, is_constant
def load_model(case_name):
pickle_path = Path.cwd().absolute().joinpath(f"{case_name}.pkl")
print(pickle_path)
assert pickle_path.exists()
assert pickle_path.is_file()
with open(pickle_path, 'rb') as f:
m = pickle.load(f)
return m
def pw_trans(m, num_points=3):
trans = pe.TransformationFactory("contrib.piecewise.nonlinear_to_pwl")
trans.apply_to(
m,
num_points=num_points,
targets=list(m.component_data_objects(pe.Constraint, descend_into=True)),
additively_decompose=True,
)
def expand_with_sympy(expr: NumericExpression):
if is_constant(expr):
return value(expr)
object_map, sympy_expr = sympyify_expression(expr, keep_mutable_parameters=True)
new_expr = sympy2pyomo_expression(sympy_expr.expand(), object_map)
if is_constant(new_expr):
new_expr = value(new_expr)
return new_expr
def expand_constraints(m):
m.expanded_constraints = pe.ConstraintList()
constraints_to_delete = []
for c in list(m.component_data_objects(pe.Constraint, descend_into=True)):
expanded_expr = expand_with_sympy(c.body)
if c.equality:
new_constraint_expr = expanded_expr == c.lower
else:
if c.lower is not None and c.upper is not None:
new_constraint_expr = pe.inequality(c.lower, expanded_expr, c.upper)
elif c.lower is not None:
new_constraint_expr = expanded_expr >= c.lower
elif c.upper is not None:
new_constraint_expr = expanded_expr <= c.upper
else:
raise ValueError(f"Constraint {c.name} has no bounds.")
# Break up linear expressions into sums of single-term expressions
new_args = [
break_up_linear_expressions(arg) for arg in new_constraint_expr.args
]
broken_up_con = new_constraint_expr.create_node_with_local_data(new_args)
m.expanded_constraints.add(broken_up_con)
constraints_to_delete.append(c)
for c in constraints_to_delete:
if c.is_indexed():
# For indexed constraints, delete the specific index
del c.parent_component()[c.index()]
else:
# For scalar constraints, delete the entire component
m.del_component(c.parent_component())
def break_up_linear_expressions(expr: SumExpression) -> SumExpression:
if not isinstance(expr, SumExpression):
return expr
new_args = []
for arg in expr.args:
if isinstance(arg, LinearExpression):
for c, v in zip(arg.linear_coefs, arg.linear_vars):
new_args.append(c * v)
else:
new_args.append(arg)
return expr.create_node_with_local_data(new_args)
if __name__ == "__main__":
m = load_model('maxmin')
expand_constraints(m) # For every constraint, simplifies using sympy. Comment or uncomment this line for a different error
pw_trans(m, num_points=3)With Triangulation.Delaunay:
ERROR: Unable to triangulate the set of input points.
ERROR: Constructing component
'_pyomo_contrib_nonlinear_to_pwl.'_pwle_expanded_constraints[1]_0'' from
data=None failed:
QhullError: QH6154 Qhull precision error: Initial simplex is flat
(facet 1 is coplanar with the interior point)
While executing: | qhull d Qz Qc Qt Qbb Q12 Options selected for Qhull
2020.2.r 2020/08/31:
run-id 115706718 delaunay Qz-infinity-point Qcoplanar-keep
Qtriangulate Qbbound-last Q12-allow-wide _pre-merge _zero-centrum
Qinterior-keep Pgood _max-width 1 Error-roundoff 2.7e-15 _one-merge
3e-14 Visible-distance 1.6e-14 U-max-coplanar 1.6e-14 Width-outside
3.3e-14 _wide-facet 9.8e-14 _maxoutside 3.3e-14
precision problems (corrected unless 'Q0' or an error)
6 nearly singular or axis-parallel hyperplanes 6 zero divisors
during back substitute
399 zero divisors during gaussian elimination
The input to qhull appears to be less than 5 dimensional, or a computation
has overflowed.
Qhull could not construct a clearly convex simplex from points:
- p2(v6): 0.0001 0.0001 0 0 0
- p1(v5): 0.0001 0.0001 0 0 0
- p81(v4): 0.5 0.5 0 0 1
- p18(v3): 0.0001 1 0 0 0.45
- p54(v2): 1 0.0001 0 0 0.45
- p0(v1): 0.0001 0.0001 0 0 0
The center point is coplanar with a facet, or a vertex is coplanar with a
neighboring facet. The maximum round off error for computing distances is
2.7e-15. The center point, facets and distances to the center point are
as follows:
center point 0.25 0.25 0 0 0.3181
facet p1 p81 p18 p54 p0 distance= 0 facet p2 p81 p18 p54 p0 distance=
0 facet p2 p1 p18 p54 p0 distance= 0 facet p2 p1 p81 p54 p0 distance=
0 facet p2 p1 p81 p18 p0 distance= 0 facet p2 p1 p81 p18 p54 distance=
0
These points either have a maximum or minimum x-coordinate, or they
maximize the determinant for k coordinates. Trial points are first
selected from points that maximize a coordinate.
The min and max coordinates for each dimension are:
0: 0.0001 0.9999 difference= 0.9998 1: 0.0001 0.9999
difference= 0.9998 2: 0 0 difference= 0 3: 0
0 difference= 0 4: 0 0.9999 difference= 0.9999
If the input should be full dimensional, you have several options that may
determine an initial simplex:
- use 'QJ' to joggle the input and make it full dimensional
- use 'QbB' to scale the points to the unit cube
- use 'QR0' to randomly rotate the input for different maximum points
- use 'Qs' to search all points for the initial simplex
- use 'En' to specify a maximum roundoff error less than 2.7e-15.
- trace execution with 'T3' to see the determinant for each point.
If the input is lower dimensional:
- use 'QJ' to joggle the input and make it full dimensional
- use 'Qbk:0Bk:0' to delete coordinate k from the input. You should
pick the coordinate with the least range. The hull will have the
correct topology.
- determine the flat containing the points, rotate the points into a
coordinate plane, and delete the other coordinates.
- add one or more points to make the input full dimensional.