3

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

我已经通过添加像 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，优化问题将按预期解决。

4

2 回答 2

3

正如您所观察到的，对变量的硬约束可能导致不可行的解决方案。我建议您通过将变量指定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

我还删除了andlb并且不给他们可能导致不可行的硬性界限。您还有一个目标是最大化with的价值。它达到最大值：当求解器能够调整值时。从图中可以看出， 的值超出了设定值范围。控制器可能无法将其保持在该范围内。约束的软约束（基于目标）方法会惩罚偏差，但不会导致不可行的解决方案。如果需要增加违反or的惩罚，参数和可以调整。ubyuTcm.Obj(-Tc)1200ySPHISPLOWSPHIWSPLO

试验结果

您似乎有一个一阶动态模型，并且您正在尝试优化 PID 参数。如果您需要对控制器输出（执行器）的饱和度进行建模，那么,if3或相应的, ,函数可能会很有用。在动态优化课程中有更多关于目标和调整的信息。max3min3if2max2min2CV

CV 调整演示

这是您的问题的可行解决方案：

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()
于 2019-09-13T13:58:11.020 回答
-1

感谢您对我的问题的广泛而有用的回答。对此，我真的非常感激。正如您正确观察到的那样，我正在尝试针对我的简单控制问题优化调整参数。我已经用软约束执行了你的代码，它确实解决了可行性问题。我还添加了 WSPHI/LO 参数并将它们的值设置得很高，以便在约束范围内找到解决方案。不过，我喜欢有一个模型，其中控制输出（“u”）有界 [0,1]。根据您的回答，我可能必须在模型中添加“if”或“max/min”语句，以避免在“u”达到界限时出现一组不可行的方程。类似于“如果 u<0, u.dt() = 0 否则 u.dt() = Kp*e ....”。是否可以添加一个变量（类型松弛变量）以确保方程组的可行性？我还将研究动态优化课程链接中的材料，以更好地理解动态建模。再次感谢您在这个问题上引导我朝着正确的方向前进。

于 2019-09-16T14:32:49.040 回答