IPOPT 被“卡住” - 编程之家

如何解决IPOPT 被“卡住”

我正在使用与 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。