如何在Gekko中选择特定时间段积分目标函数？

Gekko本身没有直接支持带时间区间参数的integral()函数，但可以通过创建时间指示参数的方式实现仅对特定时间段的变量积分进行优化。核心思路是生成一个在目标区间内为1、其余时间为0的参数，将其与待积分项相乘后再做积分，这样积分就只会在目标时间段内累积。

实现步骤 #

• 生成时间指示参数：基于模型的时间数组，创建一个参数mask，当时间在[5,8]区间内时取值为1，否则为0。
• 构造目标函数：将mask与u**2相乘后，用m.integral()计算积分，再作为目标项加入模型。

修改后的完整代码 #

#%%Import packages
import numpy as np
from gekko import GEKKO
import matplotlib.pyplot as plt

#%% Build model

#initialize GEKKO model
m = GEKKO()

#time
m.time = np.linspace(0,10,101)

#Parameters
mass1 = m.Param(value=10)
mass2 = m.Param(value=1)
final = np.zeros(np.size(m.time))
for i in range(np.size(m.time)):
    if m.time[i] >= 6.2:
        final[i] = 1
    else:   
        final[i] = 0
final = m.Param(value=final)

# 新增：创建特定时间段的指示参数（t=5到t=8）
mask = np.zeros(np.size(m.time))
for i in range(np.size(m.time)):
    if 5 <= m.time[i] <= 8:
        mask[i] = 1
    else:
        mask[i] = 0
mask = m.Param(value=mask)

#Manipulated variable
u = m.Var(value=0)

#Variables
theta = m.Var(value=0)
q = m.Var(value=0)
#Controlled Variable
y = m.Var(value=-1)
v = m.Var(value=0)

#Equations
m.Equations([y.dt() == v,
             v.dt() == mass2/(mass1+mass2) * theta + u,
             theta.dt() == q,
             q.dt() == -theta - u])

#Objective：修改为仅对5-8秒区间的u²积分进行最小化
m.Obj(final * (y**2 + v**2 + theta**2 + q**2))
m.Obj(0.001 * m.integral(mask * u**2))  # 仅目标区间积分

#%% Tuning
#global
m.options.IMODE = 6 #control

#%% Solve
m.solve()

#%% Plot solution
plt.figure()
plt.subplot(4,1,1)
plt.plot(m.time,u.value,'r-',lw=2)
plt.ylabel('Force')
# 新增：绘制指示参数，验证区间范围
plt.plot(m.time,mask.value,'k--',lw=1)
plt.legend(['u','mask'],loc='best')

plt.subplot(4,1,2)
plt.plot(m.time,v.value,'b--',lw=2)
plt.legend(['v'],loc='best')
plt.ylabel('Velocity')
plt.subplot(4,1,3)
plt.plot(m.time,y.value,'g:',lw=2)
plt.legend(['y'],loc='best')
plt.ylabel('Position')
plt.subplot(4,1,4)
plt.plot(m.time,theta.value,'m-',lw=2)
plt.plot(m.time,q.value,'k.-',lw=2)
plt.legend([r'$\theta$','q'],loc='best')
plt.ylabel('Angle')
plt.xlabel('Time')
plt.show()

补充说明 #

如果需要更平滑的区间过渡（避免阶跃指示参数可能带来的数值问题），可以用Sigmoid函数替代阶跃函数生成mask，例如：

# 用Sigmoid生成平滑过渡的mask
t = m.time
mask = 1/(1+np.exp(-10*(t-5))) - 1/(1+np.exp(-10*(t-8)))
mask = m.Param(value=mask)

这种方式在区间边界处是平滑过渡的，适合对数值稳定性要求较高的场景。

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