Python matplotlibでグラフを作る(自動車の性能曲線)

何かグラフを作ろうと考えていて、現在では自動車・バイクのカタログには性能曲線が載っていませんが、古い車・バイクには記載されていました。

ネットで下記の性能曲線を拾ってきてグラフを作ることにしました。

出力カーブの式は、性能曲線をCadに貼り付けて座標を取ります。（馬力のみ、500rpm毎程度）
それをCalc・Excelでグラフを作成し、近似曲線を作成しその時に近似式を表示させて作りました。（Calcの古いVerでは精度が悪かったですが、最近ではExcelに近い精度が出ていると思います。）

エンジン出力曲線
馬力計算の部分を関数に変更

graph.py

# !/usr/bin/python3

import matplotlib.pyplot as plt
import numpy as np
import math

maxrpm = 10000 
maxps = 50 #馬力

rpm_scale = 1000 #目盛り幅
ps_scale = 10

def ps_(rpm):
    return -2.785e-18 * rpm ** 5 + 0.00000000000004148708 * rpm ** 4 - 0.000000000100084467 * rpm ** 3 - 0.00000050100196455452 * rpm ** 2 + 0.00634708649927301 * rpm - 7.05500381512122

x = np.linspace(maxrpm / 5, maxrpm , 10000, endpoint=True)

# ps=-2.785e-18 * x ** 5 + 0.00000000000004148708 * x ** 4 - 0.000000000100084467 * x ** 3 - 0.00000050100196455452 * x ** 2 + 0.00634708649927301 * x - 7.05500381512122

ps = ps_(x)

# torque=((-2.785e-18 * x ** 5 + 0.00000000000004148708 * x ** 4 - 0.000000000100084467 * x ** 3 - 0.00000050100196455452 * x ** 2 + 0.00634708649927301 * x - 7.05500381512122)/(x*2*math.pi)*75*60)*10

torque=(ps_(x)/(x*2*math.pi)*75*60)*10 #表示の単位を10倍にしている

font = {'family':'IPAGothic'}

plt.xlim([0,maxrpm + 1000]) #グラフの範囲
plt.ylim([0,maxps + 5])

plt.title('TZR250 エンジン性能曲線', **font)  #タイトル
plt.xlabel('回転数(rpm)', **font)
plt.ylabel('トルク(*1/10 kg・m)　、　出力(PS)', **font)

plt.xticks(np.arange(0,maxrpm + 500 + rpm_scale,rpm_scale))
plt.yticks(np.arange(0,maxps + 5 + ps_scale,ps_scale))
plt.grid(color='#A4A4A4',linestyle='--')

plt.plot(x, ps, color="blue", linewidth=1, label="power")  #実線
plt.plot(x, torque, color="red", linewidth=1, label="torque")

plt.legend(loc='upper left')  #グラフの凡例の表示位置

plt.show()

問題点
凡例を日本語にしようとフォントを指定したらエラーがでてしまいました。調査の必要あり。

数式の低回転部の精度が少し悪そうです。
それ以外は結構近い値が出ています。

速度
速度計算の一部を関数に変更

graph_1.py

# !/usr/bin/python3

import matplotlib.pyplot as plt
import numpy as np
import math

maxrpm = 10000 

g1=32/13 #変速比
g2=28/16
g3=25/19
g4=26/24
g5=25/26
g6=23/27
g1g=2.37
g2g=43/14

t_dia=0.633 #タイヤ径
t_k=0.95 #タイヤ径を修正する。空気圧・車両重量による凹み・速度による遠心力等で変わる。

def speed():
    return maxrpm/(g1g*g2g)*t_dia*t_k*math.pi/1000*60

# maxspeed1=maxrpm/(g1g*g2g*g1)*t_dia*t_k*math.pi/1000*60
# maxspeed2=maxrpm/(g1g*g2g*g2)*t_dia*t_k*math.pi/1000*60
# maxspeed3=maxrpm/(g1g*g2g*g3)*t_dia*t_k*math.pi/1000*60
# maxspeed4=maxrpm/(g1g*g2g*g4)*t_dia*t_k*math.pi/1000*60
# maxspeed5=maxrpm/(g1g*g2g*g5)*t_dia*t_k*math.pi/1000*60
# maxspeed6=maxrpm/(g1g*g2g*g6)*t_dia*t_k*math.pi/1000*60

maxspeed1 = speed() / g1
maxspeed2 = speed() / g2
maxspeed3 = speed() / g3
maxspeed4 = speed() / g4
maxspeed5 = speed() / g5
maxspeed6 = speed() / g6

font = {'family':'IPAGothic'}

x = np.arange(6)
y = np.arange(2)
X, Y = np.meshgrid(x,y)
X = np.array([[0, 0, 0, 0, 0, 0], [maxspeed1, maxspeed2, maxspeed3, maxspeed4, maxspeed5, maxspeed6]])
Y = np.array([[0, 0, 0, 0, 0, 0], [maxrpm, maxrpm, maxrpm, maxrpm, maxrpm, maxrpm]])

plt.xlim([0,190]) #グラフの範囲
plt.ylim([0,maxrpm + 1000])

plt.title('TZR250 走行性能曲線', **font)  #タイトル
plt.xlabel('車速(km/h)', **font)
plt.ylabel('回転数(rpm)', **font)

plt.xticks(np.arange(0,190+20,20))
plt.yticks(np.arange(0,11000+1000,1000))
plt.grid(color='#A4A4A4',linestyle='--')

plt.plot(X, Y) 

plt.show()

駆動力・走行抵抗
一部を関数に置換え

graph_2.py

# !/usr/bin/python3

import matplotlib.pyplot as plt
import numpy as np
import math

maxrpm = 10200
g1g = 2.37 #減速比
g2g = 43 / 14
g1 = 32 / 13
g2 = 28 / 16
g3 = 25 / 19
g4 = 26 / 24
g5 = 25 / 26
g6 = 23 / 27
Transmission_efficiency = 0.88 #伝達効率(各ギヤにより異なる？）
t_dia = 0.633
t_k = 0.95

maxkmph = 190 #km/h(表示max)
Body_weight = 159 #車体重量kg
etc_weight = 61 #搭乗者、gas等kg
Rr_k = 0.015 #走行性能曲線から決定

Air_resistance_k = 0.00365 #空気抵抗係数
projected_area = 0.5 #前面投影面積m^2

Slope_k = 1.17 #勾配抵抗の係数（不明）

def slope(angle):
    return (Body_weight + etc_weight) * math.sin(angle) * Slope_k   #勾配抵抗の計算

def Rr():
    return (Body_weight + etc_weight) * Rr_k + Air_resistance_k * projected_area * x ** 2  #空気抵抗の計算

def rpm_(kmph_,gear_):
    return (kmph_ * (g1g * g2g * gear_)/(t_dia * t_k * math.pi) * 1000 / 60)  #速度から回転数を求める

def driving_force(kmph, gear):    #駆動力の計算
    return ((-2.785e-18 * rpm_(kmph, gear) ** 5 + 0.00000000000004148708 * rpm_(kmph, gear) ** 4 - 0.000000000100084467 * rpm_(kmph, gear) ** 3 - 0.00000050100196455452 * rpm_(kmph, gear) ** 2 + 0.00634708649927301 * rpm_(kmph, gear) - 7.05500381512122)/(rpm_(kmph, gear) * 2 * math.pi) * 75 * 60) * (g1g * g2g * gear) / (t_dia * t_k / 2) * Transmission_efficiency

Rolling_resistance = (Body_weight + etc_weight) * Rr_k #転がり抵抗

x = np.linspace(0, maxkmph , maxkmph, endpoint=True)
Rr0 = Rr()                           #勾配0%走行抵抗（転がり抵抗+空気抵抗）
Rr3 = Rr() + slope(math.atan(3/100)) #勾配3%走行抵抗（転がり抵抗+空気抵抗+勾配抵抗）
Rr5 = Rr() + slope(math.atan(5/100)) #勾配5%
Rr10 = Rr() + slope(math.atan(10/100)) #勾配10%
Rr15 = Rr() + slope(math.atan(15/100)) #勾配15%
Rr20 = Rr() + slope(math.atan(20/100)) #勾配20%
Rr25 = Rr() + slope(math.atan(25/100)) #勾配25%
Rr30 = Rr() + slope(math.atan(30/100)) #勾配30%
Rr35 = Rr() + slope(math.atan(35/100)) #勾配35%
Rr40 = Rr() + slope(math.atan(40/100)) #勾配40%

maxkmph = maxrpm /(g1g * g2g) * t_dia * t_k * math.pi /1000 * 60
kmph1 = np.linspace(maxkmph/g1/5, maxkmph/g1, maxkmph/g1, endpoint=True)
kmph2 = np.linspace(maxkmph/g2/5, maxkmph/g2, maxkmph/g2, endpoint=True)
kmph3 = np.linspace(maxkmph/g3/5, maxkmph/g3, maxkmph/g3, endpoint=True)
kmph4 = np.linspace(maxkmph/g4/5, maxkmph/g4, maxkmph/g4, endpoint=True)
kmph5 = np.linspace(maxkmph/g5/5, maxkmph/g5, maxkmph/g5, endpoint=True)
kmph6 = np.linspace(maxkmph/g6/5, maxkmph/g6, maxkmph/g6, endpoint=True)

driving_force1 = driving_force(kmph1, g1)
driving_force2 = driving_force(kmph2, g2)
driving_force3 = driving_force(kmph3, g3)
driving_force4 = driving_force(kmph4, g4)
driving_force5 = driving_force(kmph5, g5)
driving_force6 = driving_force(kmph6, g6)

font = {'family':'IPAGothic'}

plt.xlim([0,190]) #グラフの範囲km/h
plt.ylim([0,200]) #kgf

plt.title('TZR250 走行性能曲線', **font)  #タイトル
plt.xlabel('車速(km/h)', **font)
plt.ylabel('駆動力・走行抵抗(kgf)', **font)

plt.xticks(np.arange(0,190+20,20))
plt.yticks(np.arange(0,200+20,20))
plt.grid(color='#A4A4A4',linestyle='--')

plt.plot(x, Rr40, color="black", linewidth=1, label="40%")
plt.plot(x, Rr35, color="brown", linewidth=1, label="35%")
plt.plot(x, Rr30, color="red", linewidth=1, label="30%")
plt.plot(x, Rr25, color="blue", linewidth=1, label="25%")
plt.plot(x, Rr20, color="green", linewidth=1, label="20%")
plt.plot(x, Rr15, color="cyan", linewidth=1, label="15%")
plt.plot(x, Rr10, color="black", linewidth=1, label="10%")
# plt.plot(x, Rr10, color="black", linewidth=1, label="10%") #不要
plt.plot(x, Rr5, color="brown", linewidth=1, label="5%")
plt.plot(x, Rr3, color="red", linewidth=1, label="3%")
plt.plot(x, Rr0, color="blue", linewidth=1, label="0%")

plt.plot(kmph1, driving_force1, color="blue", linewidth=1)  #実線
plt.plot(kmph2, driving_force2, color="blue", linewidth=1)
plt.plot(kmph3, driving_force3, color="blue", linewidth=1)
plt.plot(kmph4, driving_force4, color="blue", linewidth=1)
plt.plot(kmph5, driving_force5, color="blue", linewidth=1)
plt.plot(kmph6, driving_force6, color="blue", linewidth=1)

plt.legend(loc='upper right')  #グラフの凡例の表示位置

plt.show()

走行抵抗の計算は、
http://www.isuzu.co.jp/cv/cost/manual/knowledge_2.html
を参考にしました。（単位は古いPS、kgfを使用していますが）
勾配抵抗だけが、参考と走行性能曲線と同じになりませんでしたのでSlope_kと言う係数を入れました。（勾配の角度による抵抗以外に他の何か抵抗があるのか？。他のメーカーの曲線図もいくつか見ましたが全て係数が必要で値も違うようでした。）

参考

Python matplotlibでグラフを作る（超初心者向け）
https://qiita.com/ty21ky/items/5374dec69df08584bc35
Python matplotlibでグラフを作る（超初心者向け）-2
https://qiita.com/ty21ky/items/2b1403317e769cccabac

HONDA市販ロードレースマシン「NSF250R」性能曲線

FIMロードレース世界選手権(WGP)のMoto3クラス、およびMFJ全日本ロードレース選手権のJ-GP3クラス用（NSF250Rはレース専用車のため、一般公道での走行はできません。）
http://www.honda.co.jp/HRC/NRS/SPEC/
7500rpm以下の出力と係数等は適当に決めましたので、あくまでも想像の性能曲線図です。

変更点

maxrpm = 13000 
maxps = 50 #馬力

def ps_(rpm):  #出力の計算（7500rpm以下は適当）
#    return -9.2E-019*rpm**5+0.000000000000035485*rpm**4-0.00000000053180085097*rpm**3+0.000003660002435*rpm**2-0.00562033752019423*rpm-0.780416895323216

    return 2.83E-018*rpm**5-0.000000000000114538*rpm**4+0.00000000174894095993*rpm**3-0.0000127298485578233*rpm**2+0.0503275409834266*rpm-76.8249707875472+0.2

g1=1.875
g2=1.524
g3=1.304
g4=1.167
g5=1.077
g6=1
g1g=2.952
g2g=2.333

t_dia=0.602

maxkmph = 210 #km/h(表示max)
Body_weight = 84 #車体重量kg
etc_weight = 66 #搭乗者、gas等kg
Rr_k = 0.012 #転がり抵抗係数（適当）

Air_resistance_k = 0.0028 #空気抵抗係数（適当）
projected_area = 0.5 #前面投影面積m^2（適当）

def driving_force(kmph, gear):    #駆動力の計算（7500rpm以下は適当）
#    return ((-9.2E-019*rpm_(kmph, gear)**5+0.000000000000035485*rpm_(kmph, gear)**4-0.00000000053180085097*rpm_(kmph, gear)**3+0.000003660002435*rpm_(kmph, gear)**2-0.00562033752019423*rpm_(kmph, gear)-0.780416895323216)/(rpm_(kmph, gear) * 2 * math.pi) * 75 * 60) * (g1g * g2g * gear) / (t_dia * t_k / 2) * Transmission_efficiency

    return ((2.83E-018*rpm_(kmph, gear)**5-0.000000000000114538*rpm_(kmph, gear)**4+0.00000000174894095993*rpm_(kmph, gear)**3-0.0000127298485578233*rpm_(kmph, gear)**2+0.0503275409834266*rpm_(kmph, gear)-76.8249707875472+0.2)/(rpm_(kmph, gear) * 2 * math.pi) * 75 * 60) * (g1g * g2g * gear) / (t_dia * t_k / 2) * Transmission_efficiency

これ以外にグラフの範囲等を変更する必要があります。