我想在 GEKKO 数据的非线性回归中的给定点施加值和斜率约束,请帮助我

Posted

技术标签:

【中文标题】我想在 GEKKO 数据的非线性回归中的给定点施加值和斜率约束,请帮助我【英文标题】:I want to impose value and slope constraints at a given point in non linear regression of data with GEKKO , please help me 【发布时间】:2021-11-07 14:57:18 【问题描述】:

对于这个数据,我必须执行非线性回归,但是有一些值和斜率约束,第二个 m.equation 是该点的值约束,第三个方程是斜率约束,回归器在回归过程中应该遵循这些约束,并且评估参数

import numpy as np
import matplotlib.pyplot as plt
from gekko import GEKKO  
import sympy as sp
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() #Fixed Valve single value for all data points
a.STATUS=1
b=m.FV() #Fixed Valve single value for all data points
b.STATUS=1
c=m.FV() #Fixed Valve single value for all data points
c.STATUS=1
Cp_fit=m.CV(value=Cp)  #control variable
Cp_fit.FSTATUS=1  # Feed back staus =1 \\ we tell to use the measurements

m.Equation(Cp_fit==c*T_fit**(-2)+b*T_fit+a) # model equation y=0.1*exp(a*x)
val=11.8238767562590 
slope = 0.362994963854413
e=sp.symbols('e')
m.Equation(val-((a+b*e+c*e**-2).subs(e,70)==0) 
m.Equation(slope-(sp.diff((a+b*e+c*e**-2),e).subs(e,70)==0)
 # mmodes in gekko IMODE=2 => regeression
m.options.SOLVER=1
m.solve(disp=False) # wanna se solver output
print(a.value[0],b.value[0],c.value[0])
plt.plot(T,Cp,'bo',label='data')
plt.plot(T_fit.value,Cp_fit.value,'r',label='Regression')
plt.legend()

【问题讨论】:

【参考方案1】:

问题中缺少括号。除了这个简单的修复之外,如果结果以符号方式生成然后转换为带有str() 的字符串,则可以使用 Sympy。 Gekko 变量在 Sympy 表达式中不起作用。作为字符串的 Sympy 表达式可以在 Gekko 中使用eval()。下面是一个解决d**2==7的简单例子。

from gekko import GEKKO  
import sympy as sp

m=GEKKO()

# generate symbolic equation as string with sympy
d=sp.symbols('d')
f = d**2
eqn = str(f)

# generate gekko problem
d = m.Var(1)
m.Equation(7==eval(eqn))
m.solve()
print(d.value)

您的代码类似,只需要在 Gekko 表达式之前使用 Sympy 表达式。

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

# symbolic sympy expressions
a,b,c,e=sp.symbols(('a','b','c','e'))
eqn1 = str((a+b*e+c*e**-2).subs(e,70))
eqn2 = str(sp.diff((a+b*e+c*e**-2),e).subs(e,70))

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() #Fixed Valve single value for all data points
a.STATUS=1
b=m.FV() #Fixed Valve single value for all data points
b.STATUS=1
c=m.FV() #Fixed Valve single value for all data points
c.STATUS=1
Cp_fit=m.CV(value=Cp)  #control variable
Cp_fit.FSTATUS=1  # Feed back staus =1 \\ we tell to use the measurements

m.Equation(Cp_fit==c*T_fit**(-2)+b*T_fit+a) # model equation y=0.1*exp(a*x)
val=11.8238767562590 
slope = 0.362994963854413
m.Equation(val==eval(eqn1))
m.Equation(slope==eval(eqn2))
 # mmodes in gekko IMODE=2 => regeression
m.options.SOLVER=1
m.solve(disp=True) # wanna se solver output
print(a.value[0],b.value[0],c.value[0])
plt.plot(T,Cp,'bo',label='data')
plt.plot(T_fit.value,Cp_fit.value,'r',label='Regression')
plt.legend()

这给出了求解器输出的成功解决方案:

 APMonitor, Version 1.0.1
 APMonitor Optimization Suite
 ----------------------------------------------------------------
 
 
 --------- APM Model Size ------------
 Each time step contains
   Objects      :            0
   Constants    :            0
   Variables    :            5
   Intermediates:            0
   Connections  :            0
   Equations    :            3
   Residuals    :            3
 
 Number of state variables:            593
 Number of total equations: -          826
 Number of slack variables: -            0
 ---------------------------------------
 Degrees of freedom       :           -233
 
 * Warning: DOF <= 0
 ----------------------------------------------
 Model Parameter Estimation with APOPT Solver
 ----------------------------------------------
 
 Iter    Objective  Convergence
    0  4.17639E+05  3.62995E-01
    1  9.38101E+02  2.22045E-16
    2  9.38029E+02  3.90625E-06
    3  9.38029E+02  8.32667E-17
    4  9.38029E+02  1.11022E-16
    5  9.38029E+02  1.11022E-16
    6  9.38029E+02  1.11022E-16
 Successful solution
 
 ---------------------------------------------------
 Solver         :  APOPT (v1.0)
 Solution time  :   0.171900000001187      sec
 Objective      :    938.028508697437     
 Successful solution
 ---------------------------------------------------
 
24.049369235 0.0045650595792 -61470.728583

【讨论】:

以上是关于我想在 GEKKO 数据的非线性回归中的给定点施加值和斜率约束,请帮助我的主要内容,如果未能解决你的问题,请参考以下文章

多索引 df 中的线性回归

matlab多元线性回归

第四章 线性判据与回归

SKLearn 线性回归,但在开始之前设置某些系数

吴恩达机器学习学习笔记——2.1单变量线性回归算法

个人机器学习总结之线性回归