python - 我想在 GEKKO 数据的非线性回归中的给定点施加值和斜率约束，请帮助我
问题描述
对于这个数据，我必须执行非线性回归，但有一些值和斜率约束，第二个 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()
解决方案
问题中缺少括号。除了这个简单的修复之外，如果结果以符号方式生成，然后将其转换为带有str(). 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
推荐阅读
- mongodb - 如何使用 mongo 聚合对两个或多个字段进行分组以查找差异
- ionic4 - 设置离子选择文本
- sql - Big Query - 从临时表创建表/视图
- core-data - 重命名关系时的核心数据和 CloudKit 集成问题（代码 134110）
- python-3.x - 如何在下面的 HTML 代码中使用 xpath 识别具有相同文本的多个按钮元素？
- javascript - ExtJS 7.1 Tag 字段正在执行 XSS 标签
- javascript - 给定'n'个小数位的数字的最大精度是多少
- pandas - 在圆圈内动画文本 - Matplotlib
- html - 通过 slde 将同一侧并排显示为两行
- marklogic - 如何使用 REST API 查找已安装的 marklogic 9 数据中心版本