如何用GEKKO实现带指定点值和斜率约束的非线性回归？

问题修正说明 #

你代码的核心问题是混用了SymPy符号体系和GEKKO的求解体系，GEKKO无法识别SymPy的符号运算结果，约束需要直接用GEKKO变量手动写表达式：

• 原模型为 Cp = a + b*T + c*T^(-2)，在T=70处的值约束直接代入T=70写等式即可
• 对T手动求导得斜率公式为 dCp/dT = b - 2*c*T^(-3)，代入T=70即可写斜率约束

修正后完整代码 #

import numpy as np
import matplotlib.pyplot as plt
from gekko import GEKKO  

T=np.array([  70.,   80.,   90.,  100.,  110.,  120.,  130.,  140.,  150.,
        160.,  170.,  180.,  190.,  200.,  210.,  220.,  230.,  240.,
        250.,  260.,  270.,  280.,  290.,  298.,  300.,  310.,  320.,
        330.,  340.,  343.,  350.,  360.,  363.,  370.,  380.,  383.,
        390.,  400.,  403.,  410.,  420.,  423.,  430.,  440.,  443.,
        450.,  460.,  463.,  470.,  480.,  483.,  490.,  500.,  503.,
        510.,  520.,  523.,  530.,  540.,  543.,  550.,  560.,  563.,
        570.,  580.,  583.,  590.,  600.,  610.,  620.,  623.,  630.,
        640.,  643.,  650.,  660.,  663.,  670.,  680.,  683.,  690.,
        700.,  703.,  710.,  720.,  723.,  730.,  740.,  743.,  750.,
        760.,  763.,  770.,  780.,  790.,  800.,  810.,  820.,  830.,
        840.,  850.,  860.,  870.,  880.,  890.,  900.,  910.,  920.,
        930.,  940.,  950.,  960.,  970.,  980.,  990., 1000., 1500.,
       1500.])
Cp=np.array([11.28642 , 13.19342 , 14.82796 , 16.606885, 17.3842  , 18.3733  ,
       19.21185 , 19.9262  , 20.53826 , 21.06597 , 21.52387 , 21.9238  ,
       22.27536 , 22.58634 , 22.8631  , 23.11088 , 23.33401 , 23.53603 ,
       23.71991 , 23.88818 , 24.04287 , 24.18579 , 24.31843 , 24.4     ,
       24.44204 , 24.55777 , 24.66653 , 24.7691  , 24.86624 , 24.81    ,
       24.95854 , 25.04652 , 25.02    , 25.13065 , 25.2114  , 25.24    ,
       25.28911 , 25.36401 , 25.33    , 25.43645 , 25.50675 , 25.49    ,
       25.57505 , 25.64156 , 25.6     , 25.70655 , 25.77003 , 25.7     ,
       25.83227 , 25.89344 , 25.81    , 25.95348 , 26.01259 , 26.145   ,
       26.07098 , 26.12865 , 25.98    , 26.18561 , 26.24207 , 26.04    ,
       26.29805 , 26.35354 , 26.17    , 26.4087  , 26.46352 , 26.27    ,
       26.5182  , 26.57262 , 26.62678 , 26.68089 , 26.49    , 26.73492 ,
       26.7889  , 26.59    , 26.84285 , 26.89681 , 26.69    , 26.95088 ,
       27.005   , 26.81    , 27.05915 , 27.11354 , 26.96    , 27.16812 ,
       27.22276 , 27.13    , 27.27771 , 27.33283 , 27.47    , 27.38814 ,
       27.44385 , 27.76    , 27.49973 , 27.55588 , 27.6125  , 27.66953 ,
       27.72683 , 27.78436 , 27.84238 , 27.9009  , 27.95975 , 28.01896 ,
       28.07876 , 28.13917 , 28.19976 , 28.26095 , 28.32291 , 28.38519 ,
       28.44783 , 28.51116 , 28.57536 , 28.63981 , 28.70504 , 28.77107 ,
       28.8372  , 28.90433 , 33.47658 , 33.47658 ])

m=GEKKO()
m.options.IMODE=2 
T_fit=m.Param(value=T)
a=m.FV() 
a.STATUS=1
b=m.FV() 
b.STATUS=1
c=m.FV() 
c.STATUS=1
Cp_fit=m.CV(value=Cp) 
Cp_fit.FSTATUS=1  

# 主拟合方程
m.Equation(Cp_fit == c*T_fit**(-2) + b*T_fit + a) 

val=11.8238767562590 
slope = 0.362994963854413
# T=70处的值约束
m.Equation(val == a + b*70 + c*(70**-2))
# T=70处的斜率约束，手动求导得dCp/dT = b - 2*c*T^-3
m.Equation(slope == b - 2*c*(70**-3))

m.options.SOLVER=1
m.solve(disp=False) 
print(f"拟合参数：a={a.value[0]}, b={b.value[0]}, c={c.value[0]}")

# 验证约束满足情况
pred_70 = a.value[0] + b.value[0]*70 + c.value[0]*(70**-2)
pred_slope_70 = b.value[0] - 2*c.value[0]*(70**-3)
print(f"T=70处预测值：{pred_70}，目标值：{val}，误差：{abs(pred_70-val)}")
print(f"T=70处预测斜率：{pred_slope_70}，目标斜率：{slope}，误差：{abs(pred_slope_70-slope)}")

plt.plot(T,Cp,'bo',label='原始数据')
plt.plot(T_fit.value,Cp_fit.value,'r-',label='约束回归拟合曲线')
plt.xlabel('T')
plt.ylabel('Cp')
plt.legend()
plt.show()

运行后可看到T=70处的值和斜率约束误差接近0，完全满足拟合要求。

内容的提问来源于stack exchange，提问作者Karthikeyan M mm20m013