参考链接:https://www.youtube.com/watch?v=sGZbQgDOfi4&list=PLLBUgWXdTBDjcqDl2e5F_hcBjEc6vjr1X&index=19
对已有影响石油价格的数据进行非线性建模,非线性模型为:
o
i
l
p
r
i
c
e
=
a
∗
W
T
I
b
∗
H
H
c
∗
N
G
L
d
oil_{price}=a*WTI^{b}*HH^{c}*NGL^{d}
oilprice=a∗WTIb∗HHc∗NGLd
优化目标函数为:
∑
(
o
i
l
p
r
e
d
−
o
i
l
m
e
a
s
u
r
e
d
o
i
l
m
e
a
s
u
r
e
d
)
2
\sum(\frac{oil_{pred}-oil_{measured}}{oil_{measured}})^2
∑(oilmeasuredoilpred−oilmeasured)2
数据示例如下:
利用Gekko库的模型求解程序如下:
import numpy as np
from gekko import GEKKO
import pandas as pd
import matplotlib.pyplot as plt
# 1. read txt file using pandas
# txt_file =r'http://apmonitor.com/me575/uploads/Main/oil_data.txt'
txt_file = r"D:\Files\python\gekko\oil_data.txt"
data = pd.read_csv(txt_file)
xml1 = np.array(data['WTI_PRICE'])
xml2 = np.array(data['HH_PRICE'])
xml3 = np.array(data['NGL_PRICE'])
ym = np.array(data['BEST_PRICE'])
# 2. GEKKO model
m = GEKKO(remote=False)
# 设置非线性模型中的4个变量
a = m.FV(lb=-100,ub=100)
b = m.FV(lb=-100,ub=100)
c = m.FV(lb=-100,ub=100)
d = m.FV(lb=-100,ub=100)
a.STATUS = 1
b.STATUS = 1
c.STATUS = 1
d.STATUS = 1
# 设置影响石油价格的三个变量
x1 = m.Param(value=xml1)
x2 = m.Param(value=xml2)
x3 = m.Param(value=xml3)
z = m.Param(value=ym)
y = m.Var()
# 设置约束和目标函数
m.Equation(y==a*(x1**b)*(x2**c)*(x3**d))
# m.Equation(y==a*x1**2+b*x2**2+c*x3**2+d)
m.Obj(((y-z)/z)**2)
m.options.IMODE = 2
m.options.SOLVER = 1
# 模型求解
m.solve()
print('a: ', a.value[0])
print('b: ', b.value[0])
print('c: ', c.value[0])
print('d: ', d.value[0])
cFormular = 'Formular is: '+'\n'+\
r'$A * WTI^B * HH^C * PROPANE^D$'
from scipy import stats
slope, intercept, r_value, p_value, \
std_err = stats.linregress(ym, y)
r2 = r_value**2
cR2 = "R^2 = "+str(r_value**2)
plt.figure()
plt.plot([20,140],[20,140],'k-',label='Measured')
plt.plot(ym,y,'ro',label='Predicted')
plt.xlabel('Measured Outcome (YM)')
plt.ylabel('Predicted Outcome (Y)')
plt.legend(loc='best')
plt.text(25,115,'a ='+str(a.value[0]))
plt.text(25,110,'b ='+str(b.value[0]))
plt.text(25,105,'c ='+str(c.value[0]))
plt.text(25,100,'d ='+str(d.value[0]))
plt.text(25,90, cR2)
plt.text(80,40, cFormular)
plt.grid()
plt.show()
输出结果为:
可以看到预测结果y和观察结果ym非常接近,如果将模型修改**(修改m.Equation表达式)**为以下:
o
i
l
p
r
i
c
e
=
a
∗
W
T
I
2
+
b
∗
H
H
2
+
c
∗
N
G
L
2
+
d
oil_{price}=a*WTI^{2}+b*HH^{2}+c*NGL^{2}+d
oilprice=a∗WTI2+b∗HH2+c∗NGL2+d
得到以下结果:
通过两种结果的对比可以知道非线性模型可以更好的描述石油价格与三种影响因素之间的关系。