Python 做曲线拟合和求积分的方法_Python

这是一个由加油站油罐传感器测量的油罐高度数据和出油体积，根据体积和高度的倒数，用截面积来描述油罐形状，求出拟合曲线，再用标准数据，求积分来验证拟合曲线效果和误差的一个小项目。主要的就是首先要安装Anaconda python库，然后来运用这些数学工具。

100

101

102

103

104

105

106

107

108

109

110

111

112

113

114

115

116

117

118

119

120

121

122

123

124

125

126

127

128

129

130

131

132

133

134

135

136

137

138

139

140

141

142

143

144

145

146

147

148

149

150

151

152

153

154

155

156

157

158

159

160

161

162

163

164

165

166

167

168

169

170

171

172

173

174

175

176

177

178

179

180

181

182

183

184

185

186

187

188

189

190

191

192

193

194

195

196

197

198

199

200

201

202

203

204

205

206

207

208

209

210

211

212

213

214

215

216

217

218

219

220

221

222

223

224

									###最小二乘法试验###

									import numpy as np

									import pymysql

									from scipy.optimize import leastsq

									from scipy import integrate

									###绘图，看拟合效果###

									import matplotlib.pyplot as plt

									from sympy import *

									path='E:\PythonProgram\oildata.txt'

									lieh0 =[]   #初始第一列，油管高度

									liev1 =[]   #初始第二列，油枪记录的体积

									h_median =[]  # 高度相邻中位数

									h_subtract =[]   #相邻高度作差

									v_subtract =[]   #相邻体积作差

									select_h_subtr =[]   #筛选的高度作差 ########

									select_v_subtr =[]   #筛选的体积作差

									VdtH=[]      #筛选的V 和 t 的 倒数。

									def loadData(Spath,lie0,lie1):

									 with open(Spath,'r') as f0:

									   for i in f0:

									    tmp=i.split()

									    lie0.append(float(tmp[0]))

									    lie1.append(float(tmp[2]))

									 print ("source lie0",len(lie0))

									def connectMySQL():

									 db = pymysql.connect(host='10.**.**.**', user='root', passwd='***', db="zabbix", charset="utf8") # 校罐

									 cur = db.cursor()

									 try:

									  # 查询

									  cur.execute("SELECT * FROM station_snapshot limit 10 ")

									  for row in cur.fetchall():

									   # print(row)

									   id = row[0]

									   snapshot_id = row[1]

									   DateTime = row[13]

									   attr1V = row[5]

									   attr2H = row[6]

									   print("id=%d ,snapshot_id=%s,DateTime=%s,attr1V =%s, attr2H=%s ",

									     (id, snapshot_id, DateTime, attr1V, attr2H))

									 except:

									  print("Error: unable to fecth data of station_stock")

									 try:

									  cur.execute("SELECT * FROM can_stock limit 5");

									  for row in cur.fetchall():

									   # print(row)

									   stockid = row[0]

									   stationid = row[1]

									   DateTime = row[4]

									   Volume = row[5]

									   Height = row[8]

									   print("stockid=%d ,stationid=%s,DateTime=%s,Volume =%f, Height=%f ",

									     (stockid, stationid, DateTime, Volume, Height))

									 except:

									  print("Error: unable to fecth data of can_snapshot")

									 cur.close()

									 db.close()

									def formatData(h_med,h_subtr,v_subtr):

									 lh0 = lieh0[:]

									 del lh0[0]

									 print("lh0 size(): ",len(lh0))

									 lh1 =lieh0[:]

									 del lh1[len(lh1)-1]

									 print("lh1 size() : ",len(lh1))

									 lv0 =liev1[:]

									 del lv0[0]

									 #print (liev1)

									 print ("Souce liev1 size() : ",len(liev1))

									 print ("lv1 size() :",len(lv0))

									 """

									 lv1 =liev1[:]

									 del lv1[len(lv1)-1]

									 print("lv1 size(): ",len(lv1))

									 """

									 h_med[:] =[(x+y)/2 for x,y in zip(lh0,lh1)]  ###采样点(Xi,Yi)###

									 print("h_med size（） ： ", len(h_med))

									 h_subtr[:] = [(y-x) for x,y in zip(lh0,lh1)]

									 print("h_subtr size（） ： ", len(h_subtr))

									 # v_subtr[:] = [(y-x) for x,y in zip(lv0,lv1)]

									 v_subtr[:] = lv0

									 print("v_subtr size（） ： ", len(v_subtr))

									def removeBadPoint(h_med,h_sub,v_sub):

									 for val in h_sub:

									  position=h_sub.index(val)

									  if 0.01 > val > -0.01:

									   del h_sub[position]

									   del h_med[position]

									   del v_sub[position]

									 v_dt_h_ay = [(y/x) for x, y in zip(h_sub, v_sub)]

									 return v_dt_h_ay

									def selectRightPoint(h_med,h_subtr,v_dt_h_ay):

									 for val in v_dt_h_ay:

									  pos = v_dt_h_ay.index(val)

									  if val > 20 :

									   del v_dt_h_ay[pos]

									   del h_med[pos]

									   del h_subtr[pos]

									 for val in v_dt_h_ay:

									  ptr = v_dt_h_ay.index(val)

									  if val < 14:

									   del v_dt_h_ay[ptr]

									   del h_med[ptr]

									   del h_subtr[ptr]

									def writeFile(h_mp, v_dt_h):

									 s='\n'.join(str(num)[1:-1] for num in h_mp)

									 v='\n'.join(str(vdt)[1:-1] for vdt in v_dt_h)

									 open(r'h_2.txt','w').write(s)

									 open(r'v_dt_h.txt','w').write(v)

									 print("write h_median: ",len(h_mp))

									 # print("V_dt also is (y-x) : ",v_dt_h,end="\n")

									 print("Write V_dt_h : ",len(v_dt_h))

									# file=open('data.txt','w')

									# file.write(str(h_mp));

									# file.close

									def integralCalculate(coeff,xspace):

									 vCalcute =[]

									 x = Symbol('x')

									 a, b, c, d = coeff[0]

									 y = a * x ** 3 + b * x ** 2 + c * x + d

									 i=0

									 while (i< len(xspace)-1) :

									  m = integrate(y, (x, xspace[i], xspace[i+1]))

									  vCalcute.append(abs(m))

									  i=i+1

									 print("求导结果：",vCalcute)

									 print("求导长度 len(VCalcute): ",len(vCalcute))

									 return vCalcute

									 ###需要拟合的函数func及误差error###

									def func(p,x):

									 a,b,c,d=p

									 return a*x**3+b*x**2+c*x+d #指明待拟合的函数的形状，设定的拟合函数。

									def error(p,x,y):

									 return func(p,x)-y #x、y都是列表，故返回值也是个列表

									def leastSquareFitting(h_mp,v_dt_hl):

									 p0=[1,2,6,10]  #a,b,c 的假设初始值，随着迭代会越来越小

									 #print(error(p0,h_mp,v_dt_h,"cishu")) #目标是让error 不断减小

									 #s="Test the number of iteration" #试验最小二乘法函数leastsq得调用几次error函数才能找到使得均方误差之和最小的a~c

									 Para=leastsq(error,p0,args=(h_mp,v_dt_hl)) #把error函数中除了p以外的参数打包到args中

									 a,b,c,d=Para[0]   #leastsq 返回的第一个值是a,b,c 的求解结果，leastsq返回类型相当于c++ 中的tuple

									 print(" a=",a," b=",b," c=",c," d=",d)

									 plt.figure(figsize=(8,6))

									 plt.scatter(h_mp,v_dt_hl,color="red",label="Sample Point",linewidth=3) #画样本点

									 x=np.linspace(200,2200,1000)

									 y=a*x**3+b*x**2+c*x+d

									 integralCalculate(Para,h_subtract)

									 plt.plot(x,y,color="orange",label="Fitting Curve",linewidth=2) #画拟合曲线

									 # plt.plot(h_mp, v_dt_hl,color="blue", label='Origin Line',linewidth=1) #画连接线

									 plt.legend()

									 plt.show()

									def freeParameterFitting(h_mp,v_dt_hl):

									 z1 = np.polyfit(h_mp, v_dt_hl, 6) # 第一个拟合，自由度为6

									  # 生成多项式对象

									 p1 = np.poly1d(z1)

									 print("Z1:")

									 print(z1)

									 print("P1:")

									 print(p1)

									 print("\n")

									 x = np.linspace(400, 1700, 1000)

									 plt.plot(h_mp, v_dt_hl, color="blue", label='Origin Line', linewidth=1) # 画连接线

									 plt.plot(x, p1(x), 'gv--', color="black", label='Poly Fitting Line(deg=6)', linewidth=1)

									 plt.legend()

									 plt.show()

									def main():

									 loadData(path, lieh0, liev1)

									 connectMySQL() # 读取oildata数据库

									 formatData(h_median, h_subtract, v_subtract)

									 # 去除被除数为0对应的点，并得到v 和 h 求导 值的列表

									 VdtH[:] = removeBadPoint(h_median, h_subtract, v_subtract)

									 print("h_median1：", len(h_median))

									 print("VdtH1 : ", len(VdtH))

									 # 赛选数据，去除异常点

									 selectRightPoint(h_median, h_subtract, VdtH)

									 print("h_median2：", len(h_median))

									 print("h_subtract: ", len(h_subtract))

									 print("VdtH2 : ", len(VdtH))

									 h_mp = np.array(h_median)

									 v_dt_h = np.array(VdtH)

									 writeFile(h_mp, v_dt_h)

									 # 最小二乘法作图

									 leastSquareFitting(h_mp, v_dt_h)

									 # 多项式自由参数法作图

									 freeParameterFitting(h_mp, v_dt_h)

									if __name__ == '__main__':

									 main()