在GEKKO中能否为CV在时域不同阶段设置不同SPHI/SPLO限制？

核心思路 #

直接给CV.SPHI/CV.SPLO赋值数组无法生效，因为这两个参数仅接受标量。要实现时变的软约束（目标函数驱动，而非硬约束），可以通过定义时变参数+自定义惩罚项目标函数的方式实现，完全不需要迭代求解。

具体步骤 #

1. 定义时变的上下限参数：用m.Param()将你的液位上下限数组包装为GEKKO参数，让求解器识别为时变值。
2. 替换原CV的SPHI/SPLO设置：移除原有的m.TankLevel.SPHI/m.TankLevel.SPLO赋值，改为添加自定义目标函数，对超出时变上下限的部分施加惩罚。
3. 保留CV的轨迹平滑特性（可选）：如果需要维持原有的轨迹平滑特性，可以保留TAU等参数，结合自定义惩罚项使用。

修改后的完整代码 #

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

def G1_offline(timespace=100):
    # 定义时变的液位下限
    tk_lowlimit = [37]*100       
    tk_lowlimit[40:70] = [38]*30 
    # 定义时变的液位上限（示例：全程40，可按需修改）
    tk_highlimit = [40]*100

    m = GEKKO(remote=False)
    # 将上下限转为GEKKO时变参数
    tk_low = m.Param(value=tk_lowlimit)
    tk_high = m.Param(value=tk_highlimit)

    rundown_schedule = [100]*timespace  
    rundown_schedule[40:45] = [95]*5   

    m.time = np.linspace(0, timespace-1, timespace)

    # 定义MV/DV/FV
    m.Unit1_Feed = m.MV(value=25, lb=0, ub=60, name='Unit1 Feed')
    m.Unit2_Feed = m.MV(value=27, lb=0, ub=60, name='Unit2 Feed')
    m.Fuel = m.MV(value=10, lb=0, ub=100, name='Fuel')
    m.Rundown = m.MV(name='Rundown')  # DV
    m.Efficiency = m.FV(value=0.99, lb=0.95, ub=1, name='Efficiency')                 
    m.Rundown.value = rundown_schedule

    # 定义SV/CV
    m.Flare = m.SV(value=30, lb=0, ub=100, name='Flare')
    m.TankLevel = m.CV(value=25, lb=0, ub=300, name='tklevel')

    # 定义中间变量和方程
    m.Consumers = m.MV(value=30, lb=0, ub=130, name='Consumers')
    m.Product = m.Intermediate((m.Unit1_Feed + m.Unit2_Feed)*m.Efficiency, name='Product')
    m.Balance = m.Intermediate(m.Product - m.Consumers, name='Balance')

    m.Equation(m.TankLevel.dt() == m.Balance)
    m.Equation(m.Flare == m.Rundown - (m.Unit1_Feed + m.Unit2_Feed + m.Fuel))

    # 全局选项
    m.options.IMODE = 6        # 动态控制模式（同时求解）
    m.options.NODES = 2        # 配点节点数
    m.options.SOLVER = 1       # APOPT求解器（适合混合整数/非线性优化）
    m.options.CV_TYPE = 1      # 轨迹跟踪类型（1=绝对值误差，2=平方误差）
    m.options.CTRL_UNITS = 3   # 时间单位：小时
    m.options.CTRL_TIME = 1    # 每个时间步长1小时
    m.options.REQCTRLMODE = 3  # 控制模式
    m.options.RTOL = 1e-6
    m.options.OTOL = 1e-6
    m.options.CSV_WRITE = 2

    # MV/DV状态设置
    m.Unit1_Feed.STATUS = 1
    m.Unit2_Feed.STATUS = 1
    m.Fuel.STATUS = 1
    m.Consumers.STATUS = 1
    m.Rundown.STATUS = 0       # DV固定不变
    m.Efficiency.STATUS = 0
    m.Efficiency.FSTATUS = 1

    # CV设置：保留轨迹平滑特性，移除固定SPHI/SPLO
    m.TankLevel.STATUS = 1     # 启用CV控制
    m.TankLevel.FSTATUS = 1    # 允许反馈
    m.TankLevel.TAU = 12       # 轨迹时间常数
    m.TankLevel.TR_INIT = 0    # 不重新居中轨迹
    m.TankLevel.TR_OPEN = 1    # 轨迹开放形状

    # 添加时变软约束的惩罚项（替代原SPHI/SPLO+WSPHI/WSPLO）
    WSPLO = 20  # 低于下限的惩罚权重
    WSPHI = 20  # 高于上限的惩罚权重
    # 对低于下限的部分施加惩罚：max(0, 下限 - 液位) * 权重
    m.Obj(WSPLO * m.max(tk_low - m.TankLevel, 0))
    # 对高于上限的部分施加惩罚：max(0, 液位 - 上限) * 权重
    m.Obj(WSPHI * m.max(m.TankLevel - tk_high, 0))

    # 成本函数设置
    m.Consumers.COST = -40
    m.Unit1_Feed.COST = 5
    m.Unit2_Feed.COST = 4
    m.Fuel.COST = -2

    # 移动成本（避免MV突变）
    m.Consumers.DCOST = 15
    m.Unit1_Feed.DCOST = 5
    m.Unit2_Feed.DCOST = 5
    m.Fuel.DCOST = 1

    # MV最大步长限制
    m.Consumers.DMAX = 10
    m.Unit1_Feed.DMAX = 10
    m.Unit2_Feed.DMAX = 8
    m.Fuel.DMAX = 10

    m.Consumers.MV_STEP_HOR = 1
    m.Unit1_Feed.MV_STEP_HOR = 1
    m.Unit2_Feed.MV_STEP_HOR = 1
    m.Fuel.MV_STEP_HOR = 1

    # 求解
    m.solve(GUI=False)   

    # 处理结果
    with open(m.path+'//results.json') as f:
        results = json.load(f)
        results_df = pd.DataFrame(results)
    print(results_df)

    # 绘制液位曲线及约束
    fig = plt.figure(figsize=(14,6))
    plt.plot(results_df['time'], results_df['tklevel'], color='red', label='液位')
    plt.plot(results_df['time'], tk_lowlimit, color='green', linestyle='--', label='时变下限')
    plt.plot(results_df['time'], tk_highlimit, color='blue', linestyle='--', label='时变上限')
    plt.fill_between(x=results_df['time'], y1=tk_lowlimit, y2=tk_highlimit, color='gray', alpha=0.2)

    plt.xlabel('时间（小时）')
    plt.title('动态调度优化结果')
    plt.ylabel('液位')
    plt.legend(bbox_to_anchor=(0.0, 1), loc='upper left', borderaxespad=0.5)
    plt.minorticks_on()
    plt.grid(color='b', linestyle='--', linewidth=0.5, axis='y')
    plt.show()

    # 绘制变量曲线
    fig = plt.figure(figsize=(14,6))
    plt.plot(results_df['time'], results_df['unit1_feed'], color='red', label='装置1进料')
    plt.plot(results_df['time'], results_df['unit2_feed'], color='green', label='装置2进料')
    plt.plot(results_df['time'], results_df['consumers'], color='black', label='消费需求')
    plt.plot(results_df['time'], results_df['flare'], color='orange', label='放空量')
    plt.plot(results_df['time'], results_df['fuel'], color='blue', label='燃料消耗')
    plt.plot(results_df['time'], results_df['rundown'], color='purple', label='rundown调度')
    plt.xlabel('时间（小时）'), plt.ylabel('knm3/h'), plt.title('操作变量变化')  
    plt.legend(bbox_to_anchor=(0.0, 1), loc='upper left', borderaxespad=0.5)
    plt.minorticks_on()
    plt.grid(color='b', linestyle='--', linewidth=0.5, axis='y')

    return m, results_df

# 主函数调用
c1, results_df = G1_offline(100)

关键修改说明 #

1. 时变参数定义：用m.Param()包装液位上下限数组，让求解器能识别每个时间步的约束值。
2. 自定义惩罚项：使用m.max()函数构造软约束惩罚，当液位超出上下限时，惩罚项会被加入目标函数，驱动求解器尽量让液位保持在约束范围内（但不是强制硬约束，保留优化灵活性）。
3. 保留CV轨迹特性：保留了TAU等参数，确保液位变化平滑，符合动态调度的实际需求。

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