问题描述
我正在使用与 CASadi 连接的 IPOPT 来模拟 SOC 系统,但求解器总是卡住。我知道它找到了最佳解决方案,但它没有返回最佳解决方案。当它找到解决方案时,它会不断返回目标中非常小的变化,直到达到 max_iter。
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
1480 4.4540346e+02 5.28e-03 5.42e+03 -1.7 8.26e-02 - 1.00e+00 5.00e-01h 2
1481 4.4546874e+02 3.48e-03 4.10e+03 -1.7 4.87e-02 - 1.00e+00 5.00e-01h 2
1482 4.4544443e+02 2.13e-03 3.96e+03 -1.7 3.36e-02 - 1.00e+00 5.00e-01h 2
1483 4.4553992e+02 8.01e-05 9.22e+00 -1.7 3.82e-02 - 1.00e+00 1.00e+00H 1
1484 4.4549114e+02 8.55e-04 8.52e+03 -1.7 7.03e-02 - 1.00e+00 2.50e-01f 3
1485 4.4550202e+02 1.75e-03 4.58e+03 -1.7 4.03e-02 - 1.00e+00 5.00e-01h 2
1486 4.4551714e+02 1.53e-04 7.55e+00 -1.7 4.29e-02 - 1.00e+00 1.00e+00H 1
1487 4.4548019e+02 8.06e-04 7.84e+03 -1.7 8.26e-02 - 1.00e+00 2.50e-01f 3
这是相关的代码。 我将我的功率变量分成两个部分,因为我希望“再生”率是再生率的一半。
negative = (ca.sign(power) + 1) / 2
positive = (ca.sign(power) - 1) / 2
dsoc = (-power * negative) + (0.5 * power * positive)
negative = (ca.sign(power) + 1) / 2
positive = (ca.sign(power) - 1) / 2
dsoc = (-power * negative) + (power * positive)
但是使用此代码它会变得很脏
negative = (ca.sign(power) + 1) / 2
positive = (ca.sign(power) - 1) / 2
dsoc = (-power * negative) + (0.99*power * positive)
谢谢
解决方法
可以使用松弛变量对具有不可恢复损耗的能源系统的 SOC 建模,而不是使用 sign
函数。下面是一个 IPOPT 解决方案,其存储能量的调度效率为 50%。
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.ticker as mtick
from gekko import GEKKO
m = GEKKO(remote=False)
t = np.linspace(0,24,24*3+1)
m.time = t
m.options.SOLVER = 1
m.options.IMODE = 6
m.options.NODES = 3
m.options.CV_TYPE = 1
m.options.MAX_ITER = 300
p = m.FV() # production
p.STATUS = 1
s = m.Var(100,lb=0) # storage inventory
store = m.SV() # store energy rate
vy = m.SV(lb=0) # store slack variable
recover = m.SV() # recover energy rate
vx = m.SV(lb=0) # recover slack variable
eps = 0.7
d = m.MV(-20*np.sin(np.pi*t/12)+100)
m.periodic(s)
m.Equations([p + recover/eps - store >= d,p - d == vx - vy,store == p - d + vy,recover == d - p + vx,s.dt() == store - recover/eps,store * recover <= 0])
m.Minimize(p)
m.solve(disp=True)
#%% Visualize results
fig,axes = plt.subplots(4,1,sharex=True)
ax = axes[0]
ax.plot(t,store,'C3-',label='Store Rate')
ax.plot(t,recover,'C0-.',label='Recover Rate')
ax = axes[1]
ax.plot(t,d,'k-',label='Electricity Demand')
ax.plot(t,p,'C3--',label='Power Production')
ax = axes[2]
ax.plot(t,s,'C2-',label='Energy Inventory')
ax = axes[3]
ax.plot(t,vx,label='$S_1$')
ax.plot(t,vy,label='$S_2$')
ax.set_xlabel('Time (hr)')
for ax in axes:
ax.legend(bbox_to_anchor=(1.01,0.5),\
loc='center left',frameon=False)
ax.grid()
ax.set_xlim(0,24)
loc = mtick.MultipleLocator(base=6)
ax.xaxis.set_major_locator(loc)
plt.tight_layout()
plt.show()
这是关于该问题的 additional information。