python - 两相区域 coolprop 数据的 GEKKO bspline 问题

Question

我正在尝试在 gekko 的动态热模拟中使用 coolprop。我了解到，与其直接实现 coolprop 调用，还不如使用 bspline 或 cspline 对象来让 gekko 找到导数。但是，当尝试在流体的两相区域这样做时，bspline 对象不会导致正确的结果。下面是一个最小的工作示例，其中我显示了 gekko 结果与来自 scipy 的 2d 样条曲线的比较。我正在使用样条函数根据压力和比焓计算流体的温度。

import numpy as np
from CoolProp.CoolProp import PropsSI as PS
from gekko import GEKKO
from matplotlib import pyplot as plt
from scipy.interpolate import bisplrep, bisplev


def coolprop_bspline():
    # Gekko Setup
    m = GEKKO(remote=False)
    t_steps = 90
    m.time = np.linspace(0, 90, t_steps)

    # Initial values
    xi = 0.0
    p_i = 10e5
    tau_i = PS('T', 'Q', xi, 'P', p_i, 'R134a')
    h_i = PS('H', 'Q', xi, 'P', p_i, 'R134a')
    tau_low = tau_i - 30
    tau_high = tau_i + 30

    # Prepare spline data
    h_low = PS('H', 'T', tau_low, 'P', p_i, 'R134a')
    h_high = PS('H', 'T', tau_high, 'P', p_i, 'R134a')
    h_range = np.squeeze(np.linspace(h_low, h_high, 50))
    p_range = np.linspace(9.0e5, 11e5, 10)
    h_mesh, p_mesh = np.meshgrid(h_range, p_range, indexing='ij')
    tau_mesh = np.reshape(
        PS('T', 'H', np.reshape(h_mesh, (h_mesh.size,)), 'P', np.reshape(p_mesh, (p_mesh.size,)), 'R134a'),
        h_mesh.shape)

    # Variables
    h = m.Var(h_i)
    p = m.Var(p_i)
    tau = m.Var(tau_i)

    # Equations
    m.Equations(
        [
            h.dt() == 1000,
            p.dt() == 0,
        ]
    )
    m.bspline(h, p, tau, h_range, p_range, tau_mesh, data=True, kx=3, ky=3)

    # Execution
    m.options.IMODE = 4
    m.options.SOLVER = 3
    m.solve(disp=True)

    # Evaluate comparative spline from Scipy
    spl = bisplrep(h_mesh, p_mesh, tau_mesh, kx=3, ky=3)
    h_new = np.linspace(h_low, h_high, 200)
    p_new = np.linspace(9.0e5, 11e5, 200)
    h_new, p_new = np.meshgrid(h_new, p_new, indexing='ij')
    tau_interp = bisplev(h_new[:, 0], p_new[0, :], spl)

    # Plotting

    # Plot result from Gekko
    plt.figure(1)
    plt.plot(h.value,tau.value)
    plt.grid()
    plt.xlabel('Enthalpy [J/kg]')
    plt.ylabel('Temperature [K]')
    plt.title('Gekko results with p = 1e6')

    # Plot comparative spline from scipy
    fig = plt.figure(2)
    ax = fig.add_subplot(111, projection='3d')
    ax.plot_surface(h_new, p_new, tau_interp)
    ax.view_init(20, -100)
    plt.xlabel('Enthalpy [J/kg]')
    plt.ylabel('Pressure [Pa]')
    ax.set_zlabel('Temperature [K]')
    plt.title('Correct BSpline')

    # Plot original data from coolprop
    fig = plt.figure(3)
    ax = fig.add_subplot(111, projection='3d')
    ax.plot_surface(h_mesh, p_mesh, tau_mesh)
    ax.view_init(20,-100)
    plt.xlabel('Enthalpy [J/kg]')
    plt.ylabel('Pressure [Pa]')
    ax.set_zlabel('Temperature [K]')
    plt.title('Raw coolprop data')

    plt.show()


if __name__ == '__main__':
    coolprop_bspline()

下面是结果图。将 Gekko 的结果与其他两个结果进行比较时，很明显，在 335000 J/kg 左右的焓值下，温度不应突然升高到样条覆盖的最大值。

有没有人知道如何调整 bspline 对象以正确拟合数据？到目前为止，我已经尝试了不同的多项式次数(kx、ky)以及压力和焓的对数输入。

Answer 1

Gekko 可以从数据生成 bspline，也可以使用来自另一个 bspline 生成器的结和系数。这是一个 example使用来自另一个生成器的结和系数。

from gekko import GEKKO
import numpy as np
#knots and coeffs
m = GEKKO(remote=False)
tx = [ -1, -1, -1, -1,  1,  1,  1,  1]
ty = [ -1, -1, -1, -1,  1,  1,  1,  1]
c = [1.0, 0.33333333, -0.33333333, -1.0, 0.33333333, 0.11111111, -0.11111111,
 -0.33333333, -0.33333333, -0.11111111, 0.11111111, 0.33333333, -1.0, -0.33333333,
 0.33333333, 1.0]
x = m.Var(0.5,-1,1)
y = m.Var(0.5,-1,1)
z = m.Var(2)
m.bspline(x,y,z,tx,ty,c,data=False)
m.Obj(z)
m.solve()

Gekko 还可以从数据生成 bspline

from gekko import GEKKO
import numpy as np
#raw data
m = GEKKO(remote=False)
xgrid = np.linspace(-1, 1, 20)
ygrid = xgrid
z_data = x*y
x = m.Var(0.5,-1,1)
y = m.Var(0.5,-1,1)
z = m.Var(2)
m.bspline(x,y,z,xgrid,ygrid,z_data)
m.Obj(z)
m.solve()

但如果您已经有了想要实现的东西，我建议您使用自己生成的那个。