如何在Gekko中处理向量优化

Question

我正在尝试使用Gekko这个工具来优化时间（一个向量），以便根据几个限制条件生成日程安排，主要是价格和事件发生的次数（比如我们把次数限制在每周5次）。

在我调用求和函数时，出现了以下错误：TypeError: x必须是一个包含GEKKO参数、变量或表达式的Python列表。

即使我去掉了求和的部分，我也得到了以下错误：TypeError: @error: 模型表达式 *** 函数字符串的语法错误：无效元素：ame

以下变量是长度为42的向量：
sp,
base,
neg_ln,

以下变量是标量值（长度为1）：
min,
max,
co_ln,
co1,

我可以通过一次只查看42行中的1行来优化x1，但我无法弄清楚如何在一次操作中跨多行进行优化，同时限制频率（x7）。任何指导都非常感谢。

m = GEKKO(server='XXX', remote=False)
m.options.IMODE = 2
m.options.MAX_ITER = 1000

#initialize variables
x1, x2, x3, x4, x5, x6, x7 = [m.Var() for i in range(7)]

x1.value = 1
x2.value = 1
x3.value = 1
x4.value = 0
x5.value = 0
x6.value = 0
x7.value = 1

neg_ln  =  m.Intermediate(-m.log(x1/sp))

vol1  =  m.Intermediate(co1+ base + (neg_ln*co_ln))
vol2  =  m.Intermediate(co2+ base + (neg_ln*co_ln))
vol3  =  m.Intermediate(co3+ base + (neg_ln*co_ln))
vol4  =  m.Intermediate(base + (neg_ln*co_ln))

total_vol = m.Intermediate((
(m.max2(0,base*(m.exp(vol1)-1)) * x3 +
m.max2(0,base*(m.exp(vol2)-1)) * x4 +
m.max2(0,base*(m.exp(vol3)-1)) * x5 +
m.max2(0,base*(m.exp(vol4)-1)) * x6) + base) * x7)

m.Equation(x3+x4+x5+x6 == 1) 
m.Equation(x3+x4+x5+x6 >= 0)
m.Equation(m.sum(x7)<= 5)

m.Equation(x1 >= min)
m.Equation(x1 <= max)

z1 = m.Var(lb = 3, ub = 15)
m.Equation(z1 == x1)
z2 = m.Var(lb = 3, ub = 15)
m.Equation(z2 == x2)
z3 = m.Var(lb = 0,ub = 1, integer=True)
m.Equation(z3 == x3)
z4 = m.Var(lb = 0,ub = 1, integer=True)
m.Equation(z4 == x4)
z5 = m.Var(lb = 0,ub = 1, integer=True)
m.Equation(z5 == x5)
z6 = m.Var(lb = 0,ub = 1, integer=True)
m.Equation(z6 == x6)
z7 = m.Var(lb = 0,ub = 1, integer=True)
m.Equation(z7 == x7)

m.Maximize(m.sum(simu_total_volume))

try:
m.solve(disp = True)
except:
"error"

Answer 1

试着用列表推导式来定义中间值。我在这里随便填了一些随机的常量值，因为问题中没有给出这些常量：

sp = np.random.rand(42)
base = np.random.rand(42)
co_ln = 2; co1 = 5; co2=6; co3=7

不太清楚 x 的值是否应该是向量。这个问题可以通过设置 IMODE=3 来解决，这样时间的离散化是明确的；而 IMODE=6 则是通过 m.time 来定义时间的离散化。下面是一个成功解决的示例问题。

import numpy as np
from gekko import GEKKO
m = GEKKO(remote=False)
sp = np.random.rand(42)
base = np.random.rand(42)

x1,x2 = m.Array(m.Var,2,lb=3,ub=15,value=1)
x3,x4,x5,x6,x7 = m.Array(m.Var,5,lb=0,ub=1,value=1,integer=True)
# change default values and upper bound
x3.value = 1; x7.value=1; x7.upper = 5
x = [x1,x2,x3,x4,x5,x6,x7]

n = 42
co_ln = 2; co1 = 5; co2=6; co3=7
neg_ln=[m.Intermediate(-m.log(x[0]/sp[i])) for i in range(n)]
vol1  =[m.Intermediate(co1+ base[i] + (neg_ln[i]*co_ln)) for i in range(n)]
vol2  =[m.Intermediate(co2+ base[i] + (neg_ln[i]*co_ln)) for i in range(n)]
vol3  =[m.Intermediate(co3+ base[i] + (neg_ln[i]*co_ln)) for i in range(n)]
vol4  =[m.Intermediate(base[i] + (neg_ln[i]*co_ln)) for i in range(n)]

total_vol = [m.Intermediate((
  (m.max2(0,base[i]*(m.exp(vol1[i])-1)) * x[2] +
   m.max2(0,base[i]*(m.exp(vol2[i])-1)) * x[3] +
   m.max2(0,base[i]*(m.exp(vol3[i])-1)) * x[4] +
   m.max2(0,base[i]*(m.exp(vol4[i])-1)) * x[5]) + base[i]) * x[6]) for i in range(n)]

m.Equation(sum(x[2:7])==1) 
m.Equation(sum(x[2:7])>=0)
m.Equation(x7<=5)

m.Maximize(m.sum(total_vol))

m.options.SOLVER=1
m.solve(disp = True)

for i in range(7):
    print(f'x[{i+1}]: {x[i].value[0]}')

这会产生一个解决方案：

 Number of state variables:    892
 Number of total equations: -  718
 Number of slack variables: -  2
 ---------------------------------------
 Degrees of freedom       :    172

 ----------------------------------------------
 Steady State Optimization with APOPT Solver
 ----------------------------------------------
Iter:     1 I:  0 Tm:      0.32 NLPi:   61 Dpth:    0 Lvs:    3 Obj: -4.87E+02 Gap:       NaN
--Integer Solution:  -2.10E+01 Lowest Leaf:  -4.87E+02 Gap:   1.83E+00
Iter:     2 I:  0 Tm:      0.01 NLPi:    2 Dpth:    1 Lvs:    2 Obj: -2.10E+01 Gap:  1.83E+00
--Integer Solution:  -2.10E+01 Lowest Leaf:  -4.87E+02 Gap:   1.83E+00
Iter:     3 I:  0 Tm:      0.01 NLPi:    2 Dpth:    1 Lvs:    1 Obj: -2.10E+01 Gap:  1.83E+00
Iter:     4 I:  0 Tm:      0.01 NLPi:    2 Dpth:    1 Lvs:    0 Obj:  0.00E+00 Gap:  1.83E+00
 No additional trial points, returning the best integer solution
 Successful solution

 ---------------------------------------------------
 Solver         :  APOPT (v1.0)
 Solution time  :  0.36019999999999985 sec
 Objective      :  -20.951499798000548
 Successful solution
 ---------------------------------------------------


x[1]: 3.0
x[2]: 3.0
x[3]: 0.0
x[4]: 0.0
x[5]: 0.0
x[6]: 0.0
x[7]: 1.0

这里有一些建议：

• 不要在定义 z 变量时添加整数约束。应该在声明 x 变量时就定义整数约束。
• 使用 m.Array() 来高效地创建变量数组。x 的值可以定义为一个向量 m.Array(m.Var,7) 或者一个矩阵 m.Array(m.Var,(7,42))。
• 在初始化时使用 lb 和 ub 来定义边界，或者在初始化后用 x[6].upper=5 来定义边界。使用 m.Equation(x[6]<=5) 的效率较低，因为不等式约束是隐式求解的。
• 对于未来的问题，请提供一个最小且完整的示例，这样才能重现错误。问题中缺少了几个值。