添加平方根计算后Gekko代码出现'Solution Not Found'错误求助

核心问题分析 #

Gekko非线性求解器依赖光滑的符号导数，Python原生的**(1/2)平方根写法会引入非光滑点，即便被开方数非负，也可能导致求解器无法收敛。此外，分母为零、变量初始值缺失也会加剧求解失败概率。

具体解决步骤 #

• 替换平方根写法为Gekko内置函数
  将代码中所有**(1/2)替换为Gekko提供的m.sqrt()函数，它会为求解器提供连续可导的符号表达式，避免非光滑问题。示例修改：

  # 原代码
  m.Equation(K2 == (-B2-(B2**2-C2**2+A2**2)**(1/2))/(C2-A2))
  m.Equation(rr00 == (-AA2+(AA2**2-4*BB2)**(1/2))/2)
  # 修改后
  m.Equation(K2 == (-B2 - m.sqrt(B2**2 - C2**2 + A2**2)) / (C2 - A2))
  m.Equation(rr00 == (-AA2 + m.sqrt(AA2**2 - 4*BB2)) / 2)
  

• 添加分母非零约束
  针对C2-A2作为分母的情况，添加约束确保其绝对值大于极小值，避免求解过程中出现除以零：

  m.Equation(m.abs(C2 - A2) >= 1e-6)
  

• 为变量设置初始值
  给关键变量提供合理初始猜测，帮助求解器快速定位可行解：

  theta21 = m.Var(value=-math.pi + h2/2)  # 基于初始时间点的计算值
  theta23 = m.Var(value=theta20)
  rr00 = m.Var(value=rr0)
  

• 启用调试输出定位问题
  将solve(disp=False)改为disp=True，查看求解器详细日志，确认是约束冲突、数值不稳定还是其他问题导致求解失败：

  m.solve(disp=True)
  

修改后的完整代码示例 #

from gekko import GEKKO
import math
import numpy as np
import matplotlib.pyplot as plt
import pandas as pd

m = GEKKO(remote=False)
m.time = np.linspace(0.000, 0.046875, 300)

# 参数定义
h2 = 4.907/180*math.pi
beta2 = 40/360*2*math.pi
rr1 = 292
rr2 = 398
rr3 = 130
rr0 = 429.71
rr11 = 260
rr22 = 639.5
rr33 = 260
theta20 = 112.4853/360*2*math.pi
theta25 = 2*math.pi/3*m.Param(m.time)  # 修正为Gekko Param类型
theta200 = math.pi/2
theta233 = 0

# 变量及初始值
theta21 = m.Var(value=-math.pi + h2/2)
A2 = m.Var()
B2 = m.Var()
C2 = m.Var()
K2 = m.Var()
theta23 = m.Var(value=theta20)
theta211 = theta23 + math.pi

AA2 = m.Var()
BB2 = m.Var()
rr00 = m.Var(value=rr0)
rr000 = m.Var()
rr0000 = m.Var()

# 方程定义
m.Equation(theta21 == -math.pi + h2/2 - h2*(0.43989*theta25/beta2 - 0.035014*m.sin(4*math.pi*theta25/beta2)))
m.Equation(A2 == 2*rr0*rr3*m.cos(theta20) - 2*rr1*rr3*m.cos(theta21))
m.Equation(B2 == 2*rr0*rr3*m.sin(theta20) - 2*rr1*rr3*m.sin(theta21))
m.Equation(C2 == rr0**2 + rr1**2 + rr3**2 - rr2**2 - 2*rr0*rr1*(m.cos(theta20)*m.cos(theta21)+m.sin(theta20)*m.sin(theta21)))

# 修正平方根和分母约束
m.Equation(m.abs(C2 - A2) >= 1e-6)
m.Equation(K2 == (-B2 - m.sqrt(B2**2 - C2**2 + A2**2)) / (C2 - A2))
m.Equation(theta23 == m.atan(K2)*2)

m.Equation(AA2 == 2*rr33*(m.cos(theta200)*m.cos(theta233)+m.sin(theta200)*m.sin(theta233)) - 2*rr11*(m.cos(theta200)*m.cos(theta211)+m.sin(theta200)*m.sin(theta211)))
m.Equation(BB2 == rr11**2 + rr33**2 - rr22**2 - 2*rr11*rr33*(m.cos(theta211)*m.cos(theta233)+m.sin(theta211)*m.sin(theta233)))
m.Equation(rr00 == (-AA2 + m.sqrt(AA2**2 - 4*BB2)) / 2)

# 求解配置
m.options.IMODE = 6
m.solve(disp=True)

# 结果输出
df = pd.DataFrame({'time': m.time, 'rr00': rr00.value})
df.to_csv('D:/pyyhon/2축 힘 변화/test.csv', index=None)
plt.plot(m.time, rr00.value)
plt.show()

内容的提问来源于stack exchange，提问作者Ginie Kim