python - GEKKO bspline 用于两相区域中的 Coolprop 数据的问题
问题描述
我正在尝试在 gekko 的动态热模拟中使用coolprop。我了解到,必须使用 bspline 或 cspline 对象来允许 gekko 找到导数,而不是直接实现 coolprop 调用。但是,当尝试在流体的两相区域中这样做时,bspline 对象不会导致正确的结果。下面是一个最小的工作示例,我在其中显示了与 scipy 的 2d 样条曲线相比的 gekko 结果。我正在使用样条曲线计算流体的温度,具体取决于压力和比焓。
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) 以及压力和焓的对数输入。
解决方案
Gekko 可以从数据生成 bspline,也可以使用来自另一个 bspline 生成器的节点和系数。这是使用来自另一个生成器的结和系数的示例。
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()
但如果您已经有了想要实现的东西,我建议您使用您生成的那个。
推荐阅读
- entity-framework - Ef core 正在生成外键,而我没有告诉它这样做
- python - 解决一些区间的最大总长度的 O(n log n) 动态规划问题
- azure-devops - 访问另一个项目的用户故事
- octobercms - OctoberCMS BackendAuth::getUser() 返回 NULL
- kubernetes - 在同一命名空间中的 Pod 之间连接
- java - 单击图像滑块并转到 Android Studio 中的另一个活动
- asp.net-core - @onclick:preventDefault 在 Firefox 中不起作用
- sql - 计算满足特定条件的不同值
- android - 带有 ListView 的 Android RadioButton 组
- r - 如何将 LME4 输出提取到乳胶表?