GEKKO:函数字符串语法错误

Question

我无法理解从 GEKKO 模型收到的错误消息。

就上下文而言，该模型应该优化气 Spring 辅助门的气 Spring 力和尺寸参数，以最大限度地减少运算符(operator)关闭门所需的力。我的目的是计算在 0 到 90 度之间的一系列角度所需的力，然后将所有角度的力的绝对值相加并将该值最小化。气 Spring 长度也有限制，因此优化不会产生不切实际的尺寸参数。我已经将每个角度需要计算的力和其他值定义为中间变量列表。

错误似乎与这些中间变量列表之一有关，但我不知道如何解释错误消息中的“位置:2”，更不用说“函数字符串语法错误”了。我在这里搜索了其他一些答案，但没有找到关于如何使用此位置信息对模型进行故障排除的明确答案。消息试图告诉我什么？

apm gk_model0 <br><pre> ----------------------------------------------------------------
 APMonitor, Version 1.0.1
 APMonitor Optimization Suite
 ----------------------------------------------------------------

 @error: Model Expression
 *** Error in syntax of function string: Missing operator

Position: 2
 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0
  ?

Traceback (most recent call last):
  File "C:\Users\github\general-calculations\Waterjet_Door\optimize_gas_spring_closure.py", line 80, in <module>
    m.solve()
  File "C:\Users\AppData\Local\Programs\Python\Python39\lib\site-packages\gekko\gekko.py", line 2185, in solve
    raise Exception(response)
Exception:  @error: Model Expression
 *** Error in syntax of function string: Missing operator

Position: 2
 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0
  ?

这是我正在运行的代码:

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

# Create GEKKO model
m = GEKKO()

# Fixed parameters:
Ldoor = m.Const(value=45.0)         # Length of long leg of L-shaped door
cm_factor = m.Const(value=2/3)      # factor for placing the center of mass of the door
Lcm = m.Const(value=cm_factor*Ldoor)    # Distance along door to center of mass
W = m.Const(value=250.0)            # Door weight (lbs)
ns = m.Const(value=2) # Number of gas springs in the design
wdoor = m.Const(value=5.75)         # Length of short leg of door

min_angle = 1.0
max_angle = 90.0
npts = 90
theta = np.linspace(min_angle, max_angle, num=npts, endpoint=True)

# Define thetarad as a model parameter array
thetarad = m.Array(m.Param, npts, lb=min_angle, ub=max_angle)
for i,ti in enumerate(theta):
    thetarad[i].value = ti*np.pi/180.0

# Design parameters:
ds = m.Var(28.0, lb=1.0, ub=Ldoor)
S = m.Var(200.0, lb=0.0, ub=250.0)
hp = m.Var(17.5, lb=15.0, ub=20.0)
wp = m.Var(2.0, lb=-5.0, ub=5.0)
sumabsforce = m.Var()

# Derived parameters
# S:
xs = [m.Intermediate(wdoor*m.cos(thetarad[i]) + ds*m.sin(thetarad[i])) for i in range(npts)]
ys = [m.Intermediate(ds*m.cos(thetarad[i]) - wdoor*m.sin(thetarad[i])) for i in range(npts)]
Ls = m.Intermediate(m.sqrt(wdoor**2 + ds**2))

# W:
xw = [m.Intermediate(wdoor*m.cos(thetarad[i]) + Lcm*m.sin(thetarad[i])) for i in range(npts)]
yw = [m.Intermediate(Lcm*m.cos(thetarad[i]) - wdoor*m.sin(thetarad[i])) for i in range(npts)]
Lw = m.Intermediate(m.sqrt(wdoor**2 + Lcm**2))

# F:
xf = [m.Intermediate(wdoor*m.cos(thetarad[i]) + Ldoor*m.sin(thetarad[i])) for i in range(npts)]
yf = [m.Intermediate(Ldoor*m.cos(thetarad[i]) - wdoor*m.sin(thetarad[i])) for i in range(npts)]
Lf = m.Intermediate(m.sqrt(wdoor**2 + Ldoor**2))

# Gas spring length
ls = [m.Intermediate(m.sqrt((xs - wp)**2 + (hp - ys)**2)) for i in range(npts)]

# Angles
phi = [m.Intermediate(m.atan((hp-ys[i])/(xs[i]-wp))) for i in range(npts)]
beta = [m.Intermediate(m.atan(ys[i]/xs[i])) for i in range(npts)]
gamma = [m.Intermediate(m.atan(yw[i]/xw[i])) for i in range(npts)]
kappa = [m.Intermediate(m.atan(yf[i]/xf[i])) for i in range(npts)]
alpha = [m.Intermediate(beta[i] + phi[i]) for i in range(npts)]
psi = [m.Intermediate(np.pi/2 - gamma[i]) for i in range(npts)]

# Calculate the force required at each angle:
force = [m.Intermediate((W*Lw*m.sin(psi[i]) - ns*S*Ls*m.sin(alpha[i])) / Lf) for i in range(npts)]

# Spring length constraint
maxspring = m.Intermediate(ls[0].value)
minspring = m.Intermediate(100.0)
for i,lsi in enumerate(ls):
    if lsi.value > maxspring.value:
        maxspring.value = lsi.value
    if lsi.value < minspring.value:
        minspring.value = lsi.value

m.Equation(maxspring<(2*minspring))     # Reject unrealistic gas spring extension

absforce = [m.Intermediate(m.abs(force[i])) for i in range(npts)]
m.Equation(sumabsforce==m.sum(absforce))

m.Minimize(sumabsforce)

m.options.SOLVER=1
m.solve()
print(ds, S, hp, wp)

Answer 1

定位错误的一种方法是打开运行目录并检查文本模型文件gk_model0.apm。

Model
Constants
    i0 = 45.0
    i1 = 0.6666666666666666
    i2 = 0
    i3 = 250.0
    i4 = 2
    i5 = 5.75
End Constants
Parameters
    p1 = 0.017453292519943295, <= 90.0, >= 1.0
    p2 = 0.03490658503988659, <= 90.0, >= 1.0
    p3 = 0.05235987755982988, <= 90.0, >= 1.0
    p4 = 0.06981317007977318, <= 90.0, >= 1.0
    p5 = 0.08726646259971647, <= 90.0, >= 1.0
    p6 = 0.10471975511965977, <= 90.0, >= 1.0
    p7 = 0.12217304763960307, <= 90.0, >= 1.0
    p8 = 0.13962634015954636, <= 90.0, >= 1.0
    p9 = 0.15707963267948966, <= 90.0, >= 1.0
    ...
    i546=(((i0)*(cos(p89)))-((i5)*(sin(p89))))
    i547=(((i0)*(cos(p90)))-((i5)*(sin(p90))))
    i548=sqrt((((i5)^(2))+((i0)^(2))))
    i549=sqrt((((([0, 0, 0, 0, 0, 0, 0, 0, 0, ...]))^(2))))

i548 似乎是这个等式:Lf = m.Intermediate(m.sqrt(wdoor**2 + Ldoor**2))。下一个等式是传递列表而不是数字或 gekko 变量类型的等式。

问题中间方程

ls = [m.Intermediate(m.sqrt((xs - wp)**2 + (hp - ys)**2)) \
        for i in range(npts)]

修正方程

修正后的方程有xs[i]和ys[i]。

ls = [m.Intermediate(m.sqrt((xs[i] - wp)**2 + (hp - ys[i])**2)) \
        for i in range(npts)]

通过此修正，问题仍然没有解决。将 m.abs() 替换为 m.abs3() 以获得成功的解决方案。

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

# Create GEKKO model
m = GEKKO()

# Fixed parameters:
Ldoor = m.Const(value=45.0)         # Length of long leg of L-shaped door
cm_factor = m.Const(value=2/3)      # factor for placing the center of mass of the door
Lcm = m.Const(value=cm_factor*Ldoor)    # Distance along door to center of mass
W = m.Const(value=250.0)            # Door weight (lbs)
ns = m.Const(value=2) # Number of gas springs in the design
wdoor = m.Const(value=5.75)         # Length of short leg of door

min_angle = 1.0
max_angle = 90.0
npts = 90
theta = np.linspace(min_angle, max_angle, num=npts, endpoint=True)

# Define thetarad as a model parameter array
thetarad = m.Array(m.Param, npts, lb=min_angle, ub=max_angle)
for i,ti in enumerate(theta):
    thetarad[i].value = ti*np.pi/180.0

# Design parameters:
ds = m.Var(28.0, lb=1.0, ub=Ldoor)
S = m.Var(200.0, lb=0.0, ub=250.0)
hp = m.Var(17.5, lb=15.0, ub=20.0)
wp = m.Var(2.0, lb=-5.0, ub=5.0)
sumabsforce = m.Var()

# Derived parameters
# S:
xs = [m.Intermediate(wdoor*m.cos(thetarad[i]) + ds*m.sin(thetarad[i])) for i in range(npts)]
ys = [m.Intermediate(ds*m.cos(thetarad[i]) - wdoor*m.sin(thetarad[i])) for i in range(npts)]
Ls = m.Intermediate(m.sqrt(wdoor**2 + ds**2))

# W:
xw = [m.Intermediate(wdoor*m.cos(thetarad[i]) + Lcm*m.sin(thetarad[i])) for i in range(npts)]
yw = [m.Intermediate(Lcm*m.cos(thetarad[i]) - wdoor*m.sin(thetarad[i])) for i in range(npts)]
Lw = m.Intermediate(m.sqrt(wdoor**2 + Lcm**2))

# F:
xf = [m.Intermediate(wdoor*m.cos(thetarad[i]) + Ldoor*m.sin(thetarad[i])) for i in range(npts)]
yf = [m.Intermediate(Ldoor*m.cos(thetarad[i]) - wdoor*m.sin(thetarad[i])) for i in range(npts)]
Lf = m.Intermediate(m.sqrt(wdoor**2 + Ldoor**2))

# Gas spring length
ls = [m.Intermediate(m.sqrt((xs[i] - wp)**2 + (hp - ys[i])**2)) for i in range(npts)]

# Angles
phi = [m.Intermediate(m.atan((hp-ys[i])/(xs[i]-wp))) for i in range(npts)]
beta = [m.Intermediate(m.atan(ys[i]/xs[i])) for i in range(npts)]
gamma = [m.Intermediate(m.atan(yw[i]/xw[i])) for i in range(npts)]
kappa = [m.Intermediate(m.atan(yf[i]/xf[i])) for i in range(npts)]
alpha = [m.Intermediate(beta[i] + phi[i]) for i in range(npts)]
psi = [m.Intermediate(np.pi/2 - gamma[i]) for i in range(npts)]

# Calculate the force required at each angle:
force = [m.Intermediate((W*Lw*m.sin(psi[i]) - ns*S*Ls*m.sin(alpha[i])) / Lf) for i in range(npts)]

# Spring length constraint
maxspring = m.Intermediate(ls[0].value)
minspring = m.Intermediate(100.0)
for i,lsi in enumerate(ls):
    if lsi.value > maxspring.value:
        maxspring.value = lsi.value
    if lsi.value < minspring.value:
        minspring.value = lsi.value

m.Equation(maxspring<(2*minspring))     # Reject unrealistic gas spring extension

absforce = [m.Intermediate(m.abs3(force[i])) for i in range(npts)]
m.Equation(sumabsforce==m.sum(absforce))

m.Minimize(sumabsforce)

m.options.SOLVER=1
m.open_folder()
m.solve()
print(ds.value[0], S.value[0], hp.value[0], wp.value[0])

ds、S、hp、wp的解是1.0 123.75259617 20.0 1.0445275485 。删除 m.open_folder() 以不打开运行文件夹。