Python中限制拟合函数z单调递增的3D数据插值问题

现有三维数据如下：

x = array([ 0, 0.08885313,  0.05077321,  0.05077321,  0.03807991,  0.03807991,
        0.03807991,  0.02538661,  0.02538661,  0.0126933 ,  0.0126933 ,
       -0.        , -0.        , -0.0126933 , -0.0126933 , -0.02538661,
       -0.02538661, -0.02538661, -0.03807991, -0.05077321, -0.05077321,
       -0.05077321, -0.06346652, -0.07615982, -0.11423973, -0.12693304,
       -0.13962634])
y = array([ 0, -0.15231964, -0.08885313, -0.17770625, -0.08885313, -0.10154643,
       -0.12693304, -0.07615982, -0.08885313, -0.07615982, -0.08885313,
       -0.07615982, -0.10154643, -0.08885313, -0.10154643, -0.07615982,
       -0.08885313, -0.17770625, -0.10154643, -0.08885313, -0.11423973,
       -0.12693304, -0.10154643, -0.08885313, -0.17770625, -0.12693304,
       -0.11423973])
z = array([ 0, 1.21839241, 0.78673339, 1.21839241, 0.70318648, 0.82850684,
       0.96078945, 0.64748854, 0.71014872, 0.63356406, 0.71711096,
       0.65445078, 0.77977115, 0.73103545, 0.77977115, 0.68926199,
       0.73799769, 1.19750569, 0.85635581, 0.84243133, 0.96078945,
       1.02344963, 0.93294048, 0.97471393, 1.24624138, 1.20446793,
       1.22535466])

采用以下代码进行(x,y,z)空间曲面拟合：

# 找到x的最大最小值索引及对应值
i_max, = np.where(np.isclose(x, np.max(x)))
i_min, = np.where(np.isclose(x, np.min(x)))
x_max = x[i_max][0]
x_min = x[i_min][0]

# 获取x最大最小值对应的y值
y_max = y[i_max][0]
y_min = y[i_min][0]

# 计算二维直线系数的函数
def line_coeffs(points):
    x_coords, y_coords = zip(*points)
    A = vstack([x_coords, ones(len(x_coords))]).T
    # y = a*x + b
    a, b = lstsq(A, y_coords, rcond=None)[0]
    return (a, b)

# 计算XY平面内限制所有点左右边界的两条直线系数
k1_max, k2_max = line_coeffs([(x_max, y_max), (0, 0)])
k1_min, k2_min = line_coeffs([(x_min, y_min), (0, 0)])

# 定义拟合用的数学函数
def func(xy, a, b, c, d, e, f):
    x, y = xy
    return a + b*x + c*y + d*x**2 + e*y**2 + f*x*y

# 创建曲面上的网格点
x_fill = np.linspace(-0.2, 0.2, 40)
y_fill = np.linspace(-0.3, 0, 40)
X_fill, Y_fill = np.meshgrid(x_fill, y_fill)

# 计算网格点对应的拟合z值
Z_fill = func((X_fill, Y_fill), *popt)

# 绘制数据点和拟合曲面
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.scatter(x, y, z, color='blue', label='Data Points')

ax.scatter(np.where(X_fill >= (Y_fill - k2_min)/k1_min, np.where(X_fill <= (Y_fill - k2_max)/k1_max, X_fill, np.nan), np.nan), Y_fill, Z_fill, color='red', alpha=0.5, s= 2)

ax.set_xlim([-0.2, 0.2])
ax.set_ylim([-0.3, 0.2])
ax.set_zlim([0, 1.2])

ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_zlabel('z')

plt.show()

当前拟合结果中，红色拟合点在y=-0.3附近出现z值下降，不符合预期的z值持续递增要求，需修改拟合函数使曲面满足z值单调递增特性。

要实现z值随y减小（向-0.3方向）持续递增，可通过以下几种方式修改：

1. 改用天然单调的函数形式 #

替换原二次多项式为包含y的单调递增项的函数，同时保留x方向的拟合项：

def func(xy, a, b, c, d, e):
    x, y = xy
    # 利用负指数项保证y越小（越负）z越大，x的项拟合x方向的变化
    return a + b*x + c*x**2 + d*np.exp(-e*y)

只要设置e>0，exp(-e*y)会随y减小单调递增，从而保证z值的递增特性，同时x的线性/二次项可拟合x方向的数据波动。

2. 给二次多项式添加单调性约束 #

若坚持使用二次多项式，需通过约束优化强制保证y方向的单调性：

• 对y的偏导数∂z/∂y = c + 2e*y + f*x，需在整个拟合区域（y∈[-0.3,0], x∈[-0.2,0.2]）内恒小于0（y减小即Δy为负，要Δz为正则偏导数需为负）。

使用带约束的优化器替代普通curve_fit：

from scipy.optimize import minimize

# 定义平方误差损失函数
def loss(params, x, y, z):
    a,b,c,d,e,f = params
    pred = a + b*x + c*y + d*x**2 + e*y**2 + f*x*y
    return np.sum((pred - z)**2)

# 定义单调性约束：确保所有点的y偏导数≤-0.01（避免等于0）
def constraint(params):
    a,b,c,d,e,f = params
    # 合并数据点和网格点的(x,y)
    x_all = np.concatenate([x, X_fill.flatten()])
    y_all = np.concatenate([y, Y_fill.flatten()])
    dy = c + 2*e*y_all + f*x_all
    return -np.min(dy) - 0.01  # 约束此值≥0，即最小偏导数≤-0.01

# 用原拟合结果作为初始参数
initial_guess = popt
cons = ({'type': 'ineq', 'fun': lambda p: constraint(p)})

# 执行带约束的优化
result = minimize(loss, initial_guess, args=(x,y,z), constraints=cons)
popt_constrained = result.x

# 用约束后的参数计算拟合z值
Z_fill = func((X_fill, Y_fill), *popt_constrained)

3. 截断拟合范围+单调函数 #

若y=-0.3附近无真实数据，可直接将拟合网格的y范围限制在数据覆盖的区间（原数据y最小为-0.1777），再结合单调函数：

# 修改y_fill为数据实际覆盖的区间
y_fill = np.linspace(-0.18, 0, 40)
X_fill, Y_fill = np.meshgrid(x_fill, y_fill)

这种方法既贴合数据分布，又能保证z值的递增特性。

内容的提问来源于stack exchange，提问作者jim_athon