What is the correct usage of TransformedPiecewiseLinearFunctionND in Pyomo?

Question

Background

I'm trying to use a surrogate model in Pyomo. Given a set of data labeled x, y, and z, I would like to write z as an inexpensive function of x and y.

Issue

Pyomo has tools for multivariate piecewise linear functions. See here. I setup a simple example and my function is evaluating correctly. But there there doesn't seem to be any constraint getting added for Z. I would expect the model value to be equal to the interpolated value below. I'm guessing I did something wroing in setting up TransformedPiecewiseLinearFunctionND, but I couldn't find any examples in the documentation. Any insights would be appreciated.

Code

from pyomo.core import ConcreteModel, Var, Constraint
import pyomo.environ as pe
from pyomo.core.kernel.piecewise_library.transforms_nd import (
    PiecewiseLinearFunctionND,
    TransformedPiecewiseLinearFunctionND
)
from pyomo.core.kernel.variable import (
    variable,
    variable_list
)
import pyomo.core.kernel.piecewise_library.util as util
import numpy as np
from scipy.spatial import Delaunay

npts = 100

vlist = variable_list([variable(lb=-1, ub=1),
                       variable(lb=-1, ub=1)])

tri = util.generate_delaunay(vlist, num=npts)
x, y = tri.points.T
z = np.cos(x) * np.sin(y)

model = ConcreteModel()
model.X = Var(initialize=0)
model.Y = Var(initialize=0)
model.Z = Var(initialize=999)

f = PiecewiseLinearFunctionND(tri, z)
model.g = TransformedPiecewiseLinearFunctionND(
    f=f,
    input=(model.X, model.Y),
    output=model.Z
)

def x_rule(model):
    return model.X == 0.5

def y_rule(model):
    return model.Y == 0.5

model.x_const = Constraint(rule=x_rule)
model.y_const = Constraint(rule=y_rule)

solver = pe.SolverFactory('ipopt')
solver.solve(model)

z_exact = np.cos(0.5) * np.sin(0.5)
z_interp = f([0.5, 0.5])

x_model = pe.value(model.X)
y_model = pe.value(model.Y)
z_model = pe.value(model.Z)

print(f'Z Exact: {z_exact}')
print(f'Z Interpolated: {z_interp}')
print(f'Model values (X, Y, Z): {x_model}, {y_model}, {z_model}')

Output

Z Exact: 0.42073549240394825
Z Interpolated: 0.42067082611089646
Model values (X, Y, Z): 0.5, 0.5, 999

I've also tried adding a constraint for Z manually. This produces an error:

def z_rule(model):
    return model.Z == f([model.X, model.Y])

model.z_const = Constraint(rule=z_rule)

Answer 1

You are mixing modeling components between the pyomo.kernel and pyomo.environ modeling layers. This is not supported (this page has more information).

The multi-dimensional piecewise functionality is currently only available using the pyomo.kernel interface. An example of how to use it can be found here:

https://github.com/Pyomo/pyomo/blob/main/examples/kernel/piecewise_nd_functions.py