在 Gekko Optimizer 中使用 cspline 插值

Question

我正在尝试在我的 gekko 模型中使用 cspine 插值。在这个问题中，有一个发电厂、一个蒸汽轮机和一个电网。根据用于满足电网的容量，涡轮机将具有不同的效率。我试图实现一个gekko cspline，然后让模型调用它来根据功率产生为每个时间点提供效率。我无法让它工作。这在 Gekko 中可行吗？

from gekko import GEKKO
import numpy as np
import matplotlib.pyplot as plt

# Grid demand
t = np.linspace(0, 24, 24)
e_grid = 2.5*np.sin(t/24*(2*np.pi)) + 5

# Turbine Efficiency curve based on turbine capacity
pcap = np.linspace(.1, 1, 10)  # %capacity
cap = 10*pcap
turb_eff = .75*np.sin(cap/11*np.pi)

# build model
m = GEKKO(remote=True)
m.time = t

Econs = m.Param(e_grid)
Egen = m.MV(value=5, lb=0, ub=10)  # steam production

x = m.Param(value=cap)
y = m.Var()
Turb_spline = m.cspline(x, y, cap, turb_eff)

turb_out = m.Intermediate(Egen*Turb_spline)
m.Equation(Econs == turb_out)
m.Obj(Egen)

m.options.IMODE = 5
m.options.SOLVER = 3
m.solve()

plt.plot(t, Egen.value, label='gen')
plt.plot(t, Econs.value, label='cons')
plt.xlabel('time')
plt.ylabel('Energy')
plt.legend()

我可以使用np.polyfit 来实现它。然后，我能够将多项式添加到我的 gekko 模型中并运行它以正确调整效率。我使用以下代码而不是 cspline。

Ecap = m.Intermediate(Egen/cap)
m.Equation(turb_eff == p[0]*Ecap**5 + p[1]*Ecap**4 + p[2]*Ecap**3 + p[3]*Ecap**2 + p[4]*Ecap + p[5])
turb_out = m.Intermediate(Egen*turb_eff)
m.Equation(turb_out == Econs)

我仍然想知道如何使用 cspline，以便拟合多项式无法捕获的更复杂的模型。

Answer 1

如果我正确理解了这个问题，我会像下面这样编写目标函数。

请查看代码。这仍然很粗糙，因此您可能需要对其进行修改，但我希望您能有所了解。
m.cspline 函数的输出无法保存到其他变量。您可以根据需要调用 x 或 y 变量。
变量之间的时间数组在您的代码中不匹配。

from gekko import GEKKO
import numpy as np
import matplotlib.pyplot as plt

# Grid demand
t = np.linspace(0, 50, 50)
e_grid = 2.5*np.sin(t/24*(2*np.pi)) + 5

# Turbine Efficiency curve based on turbine capacity
pcap = np.linspace(.1, 1, 10)  # %capacity
cap = 10*pcap
turb_eff = 0.75*np.sin(cap/11*np.pi)

# build model
m = GEKKO(remote=False)
m.time = t

Econs = m.Param(e_grid)
Egen = m.MV(value=5, lb=0, ub=20)  # steam production
Egen.STATUS = 1

x = m.Var()
y = m.Var()
m.Equation(x==Egen)
m.cspline(x, y, cap, turb_eff, bound_x=True)

w1 = 10
w2 = 0
w3 = 0

turb_out = m.Var()
m.Equation(turb_out == Egen*y)
m.Obj(w1*(turb_out - Econs)**2 + w2*Egen**2)

m.options.IMODE = 6
m.solve()

print(np.max(Egen.value))

plt.figure(0)
plt.plot(cap,turb_eff, 'b.', label='data')
plt.plot(x.value[1:], y.value[1:], 'ro', label='Interpolate')
plt.xlabel('x')
plt.ylabel('y')
plt.legend()

plt.figure(1)
plt.subplot(2,1,1)
plt.plot(t[1:], turb_out[1:], label='Turb_out')
plt.plot(t[1:], Egen[1:], label='Egen')
plt.plot(t[1:], Econs[1:],'r--', label='Econs')
plt.xlabel('time')
plt.ylabel('Energy')
plt.legend()

plt.subplot(2,1,2)
plt.plot(t[1:], y[1:], label='Turb_eff')
plt.legend()