IPOPT被“卡住”

Question

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

但是使用这段代码，它得到了stucs

negative = (ca.sign(power) + 1) / 2
positive = (ca.sign(power) - 1) / 2
dsoc = (-power * negative) + (0.99*power * positive)

谢谢

Answer 1

具有不可恢复损失的能源系统的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, 'C2-', label='$S_1$')
ax.plot(t, vy, 'C3--', 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。