带正弦载荷的质量-弹簧系统GEKKO优化问题求助

核心错误原因 #

你遇到的ValueError: operands could not be broadcast together with shapes (1001,) (2,)，本质是时变正弦激励的时间序列数组（长度1001）与优化变量k1/k2（标量，形状(2,)）维度不匹配——GEKKO的优化变量是标量型，无法直接和长度为N的时间数组做逐元素运算，必须换用适配优化的建模方法。

稳态正弦响应建模适配方案 #

既然基于稳态响应假设，无需求解整个时域动态过程，直接用频域分析替代时域积分，从根源解决时变激励的维度冲突：

1. 建立系统频域运动方程（矩阵形式）：
   [M]ω²{X} + [K]{X} = {F}
   其中：
   • [M]为对角质量矩阵：diag([m1, m2, m3])
   • [K]为刚度矩阵，包含待优化变量k1、k2及固定参数km：

     [k1+km, -k1,    0   ]
     [-k1,   k1+k2, -k2  ]
     [0,     -k2,   k2+km]
     

   • ω是激励角频率，{X}为稳态位移幅值向量，{F}为激励幅值向量（假设激励作用在指定质量块上，比如m2）
2. 在GEKKO中直接建模频域方程，完全规避时间序列的维度问题。

目标函数与约束实现 #

目标函数（加权组合） #

设权重系数为w（满足0 ≤ w ≤ 1），目标函数定义为：
minimize(w*(X1² + X3²) + (1-w)*(F1² + F2²))
其中：

• X1、X3：外侧质量块的稳态位移幅值，对应km弹簧的载荷km*X1、km*X3
• F1 = k1*|X1 - X2|：k1弹簧的受力幅值
• F2 = k2*|X2 - X3|：k2弹簧的受力幅值

约束条件 #

• 运动方程约束：将频域矩阵方程拆分为3个标量方程，对应每个质量块的力平衡
• 弹簧受力超限：添加F1 ≤ F_max1、F2 ≤ F_max2（F_max为弹簧许用最大受力）
• 质量块不碰撞：添加|X1 - X2| ≥ d_min、|X2 - X3| ≥ d_min（d_min为质量块间最小安全间距）

修正后的GEKKO代码示例 #

from gekko import GEKKO
import numpy as np

# 初始化模型
m = GEKKO(remote=False)

# 已知系统参数
m1 = 1.0
m2 = 1.0
m3 = 1.0
km = 10.0
omega = 2*np.pi*5.0  # 激励角频率（5Hz）
F_amp = 5.0  # 激励幅值（作用在m2上）
w = 0.6  # 权重系数，w + (1-w) = 1
F_max1 = 20.0  # k1弹簧许用最大受力
F_max2 = 20.0  # k2弹簧许用最大受力
d_min = 0.01  # 质量块最小安全间距

# 定义优化变量
k1 = m.FV(lb=1.0, ub=50.0)  # k1取值范围
k2 = m.FV(lb=1.0, ub=50.0)  # k2取值范围
k1.STATUS = 1  # 标记为参与优化的变量
k2.STATUS = 1

# 定义位移幅值变量
X1 = m.Var()
X2 = m.Var()
X3 = m.Var()

# 频域运动方程约束
m.Equation( (k1 + km - m1*omega**2)*X1 - k1*X2 == 0 )
m.Equation( -k1*X1 + (k1 + k2 - m2*omega**2)*X2 - k2*X3 == F_amp )
m.Equation( -k2*X2 + (k2 + km - m3*omega**2)*X3 == 0 )

# 计算弹簧受力幅值
F1 = m.Intermediate(k1*m.abs(X1 - X2))
F2 = m.Intermediate(k2*m.abs(X2 - X3))

# 添加约束条件
m.Equation(F1 <= F_max1)
m.Equation(F2 <= F_max2)
m.Equation(m.abs(X1 - X2) >= d_min)
m.Equation(m.abs(X2 - X3) >= d_min)

# 定义目标函数
m.Obj( w*(X1**2 + X3**2) + (1-w)*(F1**2 + F2**2) )

# 求解优化问题
m.solve(disp=True)

# 输出结果
print(f"优化后k1: {k1.value[0]:.2f}")
print(f"优化后k2: {k2.value[0]:.2f}")
print(f"X1位移幅值: {X1.value[0]:.4f}")
print(f"X3位移幅值: {X3.value[0]:.4f}")
print(f"k1弹簧受力幅值: {F1.value[0]:.2f}")
print(f"k2弹簧受力幅值: {F2.value[0]:.2f}")

关键说明 #

• 频域建模既符合稳态响应假设，又彻底解决了时变激励的维度冲突问题，计算效率远高于时域仿真
• 使用GEKKO内置的m.abs()函数处理绝对值，避免非光滑性对优化求解的影响
• 优化变量k1、k2用FV类型定义，并设置STATUS=1标记为参与优化的变量

内容的提问来源于stack exchange，提问作者Stian HS