Gekko控制优化中运动方程传递报TypeError问题求助

我来帮你梳理并解决这两个问题，重点先搞定那个类型错误：

一、优化接触开关c的实现 #

你当前用循环生成c的方式虽然没报错，但可以用Gekko的向量化操作简化，还能避免后续和变量运算的类型问题。另外注意不要用abs做变量名，它是Python内置函数，会覆盖导致潜在问题：

# 替换原来的循环代码，直接生成Gekko数组型Intermediate
difference = m.Intermediate(0.5 - m.time)
abs_diff = m.if3(difference, -difference, difference)
c = m.Intermediate(m.if3(abs_diff - (0.5 - a), 1, 0))

这样生成的c是对应每个时间点的Gekko数组，完全符合你预期的[0,0,0,1,...,1,0,0,0]格式。

二、修复TypeError: 'int' object is not subscriptable报错 #

这个错误的核心原因是你用了Python列表c和Gekko变量做运算——列表是普通Python对象，不能直接和Gekko的MV/Var变量进行元素级运算。另外代码里还有几个小bug需要一起修正：

关键修正点： #

1. 把c改成Gekko数组：用上面优化后的c代替原来的Python列表。
2. 修正变量名错误：你写了m.fix_final(xdot,val=ydotf)，这里应该是ydot，不是xdot。
3. 修正物理公式错误：-g/m应该是-g/M，M才是你定义的质量参数。
4. 避免除以零：运动方程分母里的正弦项可能出现零值，加个小偏移量1e-6防止报错。
5. 优化控制输入定义：如果需要每个时间点都有独立的控制量，把u定义为数组型MV：u = m.MV(value=np.zeros(N)); u.STATUS = 1。

修正后的运动方程代码： #

# 替换原来的状态空间模型代码
m.Equation(xdot.dt()/TF == c*u*(L1*m.sin(q1)-L2*m.sin(q1+q2))/(M*L1*L2*m.sin(2*q1+q2)+1e-6))
m.Equation(ydot.dt()/TF == -c*u*(L1*m.cos(q1)+L2*m.cos(q1+q2))/(M*L1*L2*m.sin(2*q1+q2)+1e-6) - g/M)
m.Equation(x.dt()/TF == xdot)
m.Equation(y.dt()/TF == ydot)
m.Equation(y*init == y*final) # 初始和最终y位置相等的约束

# 修正后的终值约束
m.fix_final(x,val=xf)
m.fix_final(xdot,val=xdotf)
m.fix_final(ydot,val=ydotf)

额外建议： #

用m.fix_initial设置初始条件比直接赋值x.value = x0更清晰规范，比如：

m.fix_initial(x, val=x0)
m.fix_initial(xdot, val=xdot0)
m.fix_initial(ydot, val=ydot0)

完整修正后的核心片段 #

把这些修改整合到你的代码里，应该就能顺利传递运动方程并求解了：

# ... 前面的参数定义部分保持不变 ...

#State Variables
x, y, xdot, ydot, q1, q2 = m.Array(m.Var, 6)

# 用fix_initial设置初始条件
m.fix_initial(x, val=x0)
m.fix_initial(xdot, val=xdot0)
m.fix_initial(ydot, val=ydot0)
y.LOWER = ymin; y.UPPER = ymax
q1.LOWER = q1min; q1.UPPER = q1max
q2.LOWER = q2min; q2.UPPER = q2max

#Intermediates - 优化后的接触开关c
difference = m.Intermediate(0.5 - m.time)
abs_diff = m.if3(difference, -difference, difference)
c = m.Intermediate(m.if3(abs_diff - (0.5 - a), 1, 0))

# 控制输入（如果需要每个时间点独立控制，改成数组型MV）
u = m.MV(value=np.zeros(N)); u.STATUS = 1

#Defining the State Space Model
m.Equation(xdot.dt()/TF == c*u*(L1*m.sin(q1)-L2*m.sin(q1+q2))/(M*L1*L2*m.sin(2*q1+q2)+1e-6))
m.Equation(ydot.dt()/TF == -c*u*(L1*m.cos(q1)+L2*m.cos(q1+q2))/(M*L1*L2*m.sin(2*q1+q2)+1e-6) - g/M)
m.Equation(x.dt()/TF == xdot)
m.Equation(y.dt()/TF == ydot)
m.Equation(y*init == y*final)

#Defining final condition
m.fix_final(x,val=xf)
m.fix_final(xdot,val=xdotf)
m.fix_final(ydot,val=ydotf)

# ... 目标函数和求解部分保持不变 ...

内容的提问来源于stack exchange，提问作者T Bounds