基于GEKKO的状态切换检测求助：时变参数估计失败

我们来自温德斯海姆应用科学大学，在使用GEKKO Python进行每日多次变化的参数学习（状态切换检测）时遇到问题。现有三类时间序列数据：

• CO₂浓度（co2__ppm）
• 通风率（由阀门开度分数valve_frac与已知最大通风率相乘得到）
• 室内人数（occupancy__p）

研究问题 #

• 问题A：已知CO₂浓度和室内人数的时序数据，需估计每日每小时的阀门开度分数
• 问题B：已知CO₂浓度和通风率，估计室内人数
  当前聚焦问题A

提出的三种GEKKO实现策略 #

1. 策略I：将valve_frac声明为操纵变量（m.MV），利用MV可在每个时间点由优化器调整的特性
2. 策略II：将时间划分为短时段（如每小时一段），为每个时段声明一个固定变量（m.FV）来学习valve_frac
3. 策略III：将问题重构为控制问题，以特定CO₂浓度为设定值，valve_frac作为被控变量（m.CV）

时序数据示例 #

尝试的策略I代码 #

from gekko import GEKKO
import numpy as np
import pandas as pd

# 假设df_learn、duration__s、step__s、vent_max__m3_s_1、room__m3已定义
# Gekko Model - Initialize
m = GEKKO(remote = False)
m.time = np.arange(0, duration__s, step__s)

# Conversion factors
s_min_1 = 60
min_h_1 = 60
s_h_1 = s_min_1 * min_h_1
mL_m_3 = 1e3 * 1e3
million = 1e6

# Constants
MET__mL_min_1_kg_1_p_1 = 3.5                                  
desk_work__MET = 1.5                                  
P_std__Pa = 101325                                             
R__m3_Pa_K_1_mol_1 = 8.3145                                    
T_room__degC = 20.0                                            
T_std__degC = 0.0                                              
T_zero__K = 273.15                                            
T_std__K = T_zero__K + T_std__degC                             
T_room__K = T_zero__K + T_room__degC
infilt__m2 = 0.001                          

# Approximations
room__mol_m_3 = P_std__Pa / (R__m3_Pa_K_1_mol_1 * T_room__K)  
std__mol_m_3 = P_std__Pa / (R__m3_Pa_K_1_mol_1 * T_std__K)    
co2_ext__ppm = 415                                             
        
# National averages
weight__kg = 77.5                                              
MET__m3_s_1_p_1 = MET__mL_min_1_kg_1_p_1 * weight__kg / (s_min_1 * mL_m_3)
MET_mol_s_1_p_1 = MET__m3_s_1_p_1 * std__mol_m_3           
co2_o2 = 0.894                                            
co2__mol0_p_1_s_1 = co2_o2 * desk_work__MET * MET_mol_s_1_p_1    

# Room averages
wind__m_s_1 = 3.0                                              

# GEKKO Manipulated Variables: measured values
occupancy__p = m.MV(value = df_learn.occupancy__p.values)
occupancy__p.STATUS = 0; occupancy__p.FSTATUS = 1

# Strategy I:
valve_frac__0 = m.MV(value = df_learn.valve_frac__0.values)
valve_frac__0.STATUS = 1; valve_frac__0.FSTATUS = 0

# Strategy II:
#valve_frac__0 = m.FV(value = df_learn.valve_frac__0.values)
#valve_frac__0.STATUS = 1; valve_frac__0.FSTATUS = 0

# GEKKO Control Variable (predicted variable)
co2__ppm = m.CV(value = df_learn.co2__ppm.values) 
co2__ppm.STATUS = 1; co2__ppm.FSTATUS = 1

# GEKKO - Equations
co2_loss__ppm_s_1 = m.Intermediate((co2__ppm - co2_ext__ppm) * (vent_max__m3_s_1 * valve_frac__0 + wind__m_s_1 * infilt__m2) / room__m3)

co2_gain_mol0_s_1 = m.Intermediate(occupancy__p * co2__mol0_p_1_s_1 / (room__m3 * room__mol_m_3))

co2_gain__ppm_s_1 = m.Intermediate(co2_gain_mol0_s_1 * million)

m.Equation(co2__ppm.dt() == co2_gain__ppm_s_1 - co2_loss__ppm_s_1)

# GEKKO - Solver setting
m.options.IMODE = 5
m.options.EV_TYPE = 1
m.options.NODES = 2
m.solve(disp = False) 

执行结果异常 #

• 策略I：无模拟co2__ppm输出，valve_frac__0结果全为0
• 策略II：模拟与实测co2__ppm差异显著，valve_frac__0得到不合理的0.166（固定值）

咨询问题 #

1. 我们的GEKKO代码是否存在明显错误？
2. 策略I是否不可行，是否应切换至策略II或III？
3. GEKKO Python是否适合此类状态切换检测任务？

1. 代码中的明显错误 #

• 关键变量未定义：代码中使用vent_max__m3_s_1、room__m3、df_learn等变量但未给出定义，会直接导致运行错误或结果异常。
• 变量类型误用：
  • 实测数据occupancy__p应该用m.Param()而非m.MV()，MV是用于主动控制的操纵变量，不适合存储固定实测值。
  • 策略I中valve_frac__0作为MV未设置0-1的上下限，优化器可能取到无意义的极端值（如0）。
• CV配置缺失：co2__ppm作为CV未设置拟合权重（co2__ppm.WSP），优化器对拟合误差的重视程度不足，导致结果偏离实测值。
• MV平滑约束缺失：未给valve_frac__0设置DCOST（移动惩罚项），优化器可能让MV突变到极端值，不符合实际阀门操作的平滑性。

2. 策略可行性分析 #

• 策略I：并非不可行，需修正配置：
  • 将occupancy__p改为m.Param()
  • 给valve_frac__0设置LOWER=0、UPPER=1的上下限，添加DCOST=1e-4抑制频繁大幅变化
  • 给co2__ppm设置WSP=100，让拟合误差成为优化核心目标
• 策略II：可行，但需优化配置：
  • 确保时段划分与数据采样频率对齐
  • 给每个时段的FV设置0-1上下限，添加相邻时段的平滑约束（如m.Equation(valve_frac[i] - valve_frac[i-1] <= 0.1)）
  • 调整CV的拟合权重，提升模拟精度
• 策略III：不适合当前问题A，该策略是将CO₂作为设定值控制阀门，属于正向控制问题，与当前“已知CO₂估计阀门开度”的反向拟合需求不符。

3. GEKKO是否适合状态切换检测 #

GEKKO完全适合这类状态切换检测任务：

• 它支持动态优化与参数估计，能处理带约束的时序拟合问题。
• 针对状态切换场景，可通过以下方式实现：
  • 用分段FV（策略II）表示不同时段的参数，添加平滑约束模拟切换合理性
  • 用m.if3()等逻辑函数检测切换点，结合MV/FV实现分段参数学习
  • 采用参数自适应的MPC模式（IMODE=6），但当前问题更适合参数估计模式（IMODE=5）

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