Using a cspline interpolation in Gekko Optimizer

Question

I am trying to use a cspine interpolation in my gekko model. In this problem there is a power plant, a steam turbine, and a grid. The turbine will have different efficiencies based on the amount of it's capacity being used to meet the grid. I have tried to implement a gekko cspline and then have the model call that to give the efficiency for each time point based on the power production. I haven't been able to get this to work. Is this possible in 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()

I was able to get this to work using np.polyfit. I then was able to add a Polynomial into my gekko model and run it to get the correct adjusting of the efficiencies. I used the following code instead of a 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)

I would still like know how to use a cspline so that I can fit more complicated models that a polynomial is unable to capture.

Answer 1

If I understand the problem correctly, I would write the Objective function like below.

Please see the code. This is still very crude so that you might need to modify it, but I hope you can get an idea.
The output of m.cspline function cannot be saved to the other variable. You can just call the x or y variable as you need.
The time arrays between variables don't match in your code.

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()