gekko - 如何在 GEKKO 中约束模型变量

Question

我喜欢在 y 模型中约束变量值 u < 1。将 ub=1 添加到变量定义 u = m.Var(name='u', value=0, lb=-2, ub=1) 但它导致“找不到解决方案”(退出:收敛到一个点local infeasibility. 问题可能是不可行的。)我想我必须重新制定问题以避免这种情况，但我无法找到应该如何完成的示例。在约束变量值时，如何编写合适的模型来避免不可行的解决方案？

我已经通过添加像 m.Equation(u < 1) 这样的方程来重新表述问题，但没有成功。

import numpy as np
from gekko import GEKKO
import matplotlib.pyplot as pyplt


m = GEKKO(remote=False)

t = np.linspace(0, 1000, 101)  # time 
d = np.ones(t.shape)
d[0:10] = 0
y_delay=0

# Add data to model
m.time = t
K = m.Const(0.01, name='K')
r = m.Const(name='r', value=0)  # Reference
d = m.Param(name='d', value=d)  # Disturbance
y = m.Var(name='y', value=0, lb=-2, ub=2)  # State variable
u = m.Var(name='u', value=0, lb=-2, ub=1)  # Output
e = m.Var(name='e', value=0)
Tc = m.FV(name='Tc', value=1200, lb=60, ub=1200)  # time constant

# Update variable status
Tc.STATUS = 1  # Optimizer can adjust value

Kp = m.Intermediate(1 / K * 1 / Tc, name='Kp')
Ti = m.Intermediate(4 * Tc, name='Ti')

# Model equations
m.Equations([y.dt() == K * (u-d),
             e == r-y,
             u.dt() == Kp*e.dt()+Kp/Ti*e])

# Model constraints
m.Equation(y < 0.5)
m.Equation(y > -0.5)

# Model objective
m.Obj(-Tc)

# options
m.options.IMODE = 6  # Problem type: 6 = Dynamic optimization

# solve
m.solve(disp=True, debug=True)
print('Tc: %6.2f [s]' % (Tc.value[-1], ))

fig1, (ax1, ax2, ax3) = pyplt.subplots(3, sharex='all')
ax1.plot(t, y.value)
ax1.set_ylabel("y", fontsize=8), ax1.grid(True, which='both')
ax2.plot(t, e.value)
ax2.set_ylabel("e", fontsize=8), ax2.grid(True, which='both')
ax3.plot(t, u.value)
ax3.plot(t, d.value)
ax3.set_ylabel("u and d", fontsize=8), ax3.grid(True, which='both')
pyplt.show()

退出:收敛到局部不可行点。问题可能是不可行的。

发生错误。 错误码为2

如果我将 u 的上限更改为 2，优化问题将按预期解决。

Answer 1

如您所见，对变量的严格约束可能导致解决方案不可行。我建议您通过将变量 y 指定为受控变量并使用 SPHI 和 SPLO 设置上限和下限设置点范围来使用软约束.

y = m.CV(name='y', value=0)  # Controlled variable
y.STATUS = 1
y.TR_INIT = 0
y.SPHI = 0.5
y.SPLO = -0.5

我还从 y 和 u 中删除了 lb 和 ub导致不可行。您还有一个目标是使用 m.Obj(-Tc) 最大化 Tc 的值。它达到最大限制:1200 当求解器能够调整该值时。从图中可以看出，y 的值超出了设定点范围。 Controller 可能无法将其保持在该范围内。软约束(基于目标)的约束方法会惩罚偏差，但不会导致不可行的解决方案。如果您需要增加对违反SPHI 或SPLO 的处罚，参数WSPHI 和WSPLO 可以是调整。

看来您有一个一阶动态模型，并且您正在尝试优化 PID 参数。如果您需要模拟 Controller 输出(执行器)的饱和度，则 if3、max3、min3 或相应的 if2、max2、min2 函数可能会有用。 CV objectives 有更多信息和 tuning在动态优化类(class)中。

以下是解决您的问题的可行方案:

import numpy as np
from gekko import GEKKO
import matplotlib.pyplot as pyplt

m = GEKKO() # remote=False

t = np.linspace(0, 1000, 101)  # time 
d = np.ones(t.shape)
d[0:10] = 0
y_delay=0

# Add data to model
m.time = t
K = m.Const(0.01, name='K')
r = m.Const(name='r', value=0)  # Reference
d = m.Param(name='d', value=d)  # Disturbance
e = m.Var(name='e', value=0)
u = m.Var(name='u', value=0)  # Output
Tc = m.FV(name='Tc', value=1200, lb=60, ub=1200)  # time constant

y = m.CV(name='y', value=0)  # Controlled variable
y.STATUS = 1
y.TR_INIT = 0
y.SPHI = 0.5
y.SPLO = -0.5

# Update variable status
Tc.STATUS = 1  # Optimizer can adjust value

Kp = m.Intermediate((1 / K) * (1 / Tc), name='Kp')
Ti = m.Intermediate(4 * Tc, name='Ti')

# Model equations
m.Equations([y.dt() == K * (u-d),
             e == r-y,
             u.dt() == Kp*e.dt()+(Kp/Ti)*e])

# Model constraints
#m.Equation(y < 0.5)
#m.Equation(y > -0.5)

# Model objective
m.Obj(-Tc)

# options
m.options.IMODE = 6  # Problem type: 6 = Dynamic optimization
m.options.SOLVER = 3
m.options.MAX_ITER = 1000

# solve
m.solve(disp=True, debug=True)
print('Tc: %6.2f [s]' % (Tc.value[-1], ))

fig1, (ax1, ax2, ax3) = pyplt.subplots(3, sharex='all')
ax1.plot(t, y.value)
ax1.plot([min(t),max(t)],[0.5,0.5],'k--')
ax1.plot([min(t),max(t)],[-0.5,-0.5],'k--')
ax1.set_ylabel("y", fontsize=8), ax1.grid(True, which='both')
ax2.plot(t, e.value)
ax2.set_ylabel("e", fontsize=8), ax2.grid(True, which='both')
ax3.plot(t, u.value)
ax3.plot(t, d.value)
ax3.set_ylabel("u and d", fontsize=8), ax3.grid(True, which='both')
pyplt.show()