在 Gekko Optimizer

问题描述

我正在尝试在我的壁虎模型中使用 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 使其工作。然后，我能够将多项式添加到我的壁虎模型中并运行它以正确调整效率。我使用了以下代码而不是 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 来拟合多项式无法捕获的更复杂的模型。

解决方法

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

请看代码。这还是很粗糙的，可能需要修改一下，但希望你能有想法。
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,'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)
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,2)
plt.plot(t[1:],y[1:],label='Turb_eff')
plt.legend()