Multivariate piecewise

I don't think the current piecewise functionality supports multivariate piecewise (f(x, y)), as described in e.g., this paper.

Pyomo handles this by requesting the user to first discretise the domain using Delaunay triangulation (see here) and then to provide the y values at each of those mesh points. An example implementation is in this forum post:

import pyomo.environ as pyo
import pyomo.kernel as pmo
import numpy as np
from pyomo.core.kernel.piecewise_library.transforms_nd import piecewise_nd

model = pmo.block()
model.x = pmo.variable(lb=0, ub=3)
model.y = pmo.variable(lb=0, ub=3)
model.z = pmo.variable()

def f(x, y):
    return x*y**2 + 20*np.sin(x*y)

var_list = [model.x, model.y]
tri = pmo.piecewise_util.generate_delaunay(var_list, num=20)
values = f(*tri.points.T)

model.pw = piecewise_nd(tri, values, input=var_list, output=model.z)
model.obj = pmo.objective(expr=model.z)

pyo.SolverFactory('glpk').solve(model, tee=True)

Implementation

Assuming we limit this to a 3D space, this would require the piecewise functionality to accept an additional set of breakpoints and an additional decision variable, and the y_breakpoints to be 2D.
You could have x_points and x1_points (or z_points?).

Why

An example of when this would be useful is for operating costs in an energy system model, where you want to be able to define a cost as a function of load rate (dispatch / nominal capacity) whilst nominal capacity is a decision variable. So, you would have load rate as x, nominal capacity as z, and operating cost as y.

No source file or test is named. Start by locating the current piecewise functionality and its breakpoint and decision-variable interface, then compare it with Pyomo's piecewise_nd and the linked forum example. Done should mean supporting multivariate inputs with multidimensional breakpoints and an additional decision variable for the output.