這是一個由加油站油罐傳感器測量的油罐高度數據和出油體積，根據體積和高度的倒數，用截面積來描述油罐形狀，求出擬合曲線，再用標準數據，求積分來驗證擬合曲線效果和誤差的一個小項目。 主要的就是首先要安裝Anaconda  python庫，然后來運用這些數學工具。

###最小二乘法試驗###import numpy as npimport pymysqlfrom scipy.optimize import leastsqfrom scipy import integrate###繪圖，看擬合效果###import matplotlib.pyplot as pltfrom 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()