加工条件をデータで決める：非線形最適化で品質・電力・制約を両立する10本ノック

架空の樹脂成形工程を題材に、品質損失と電力費を抑えながら設備制約を守る条件設定を扱います。対象は No.041〜No.050（非線形最適化） です。

[!NOTE] 本資料は、数理工房 (もしくは代表である和山個人) が過去に企業研修において使用した notebook を企業様の許可を得て再構成・編集のうえ公開しています。
掲載データはすべて架空のものであり、実在する企業・工場・数値とは一切関係ありません。

はじめに：この記事で扱う製造業の実務課題

成形温度と射出圧力を決め、品質、電力、設備制約を両立します。最適値だけでなく収束、感度、制約余裕も確認します。

現場でよくある状況

試作回数には限度があり、1因子ずつの調整では曲率や相互作用を見落とします。

なぜこの問題は判断が難しいのか

改善方向が現在地で変わり、複数KPIと制約が競合するためです。数学的な成功表示と現場で安全に使えることは同義ではありません。

今回扱うノックの全体像

No.	テーマ	意思決定上の視点
041	非線形最適化	目的・変数・制約への翻訳
042	勾配	改善方向と局所感度
043	勾配降下法	歩幅と収束
044	ニュートン法	曲率の利用
045	準ニュートン法	曲率の近似
046	L-BFGS	大規模・省メモリ
047	ラグランジュ法	等式制約
048	KKT条件	制約付き解の点検
049	内点法	実行可能領域内の探索
050	機械学習との関係	予測から条件推薦へ

Python 環境の準備

外部データは使わず、乱数シードを固定します。グラフはMatplotlibで作成します。

import platform, numpy as np, pandas as pd, matplotlib, matplotlib.pyplot as plt, scipy
from scipy.optimize import minimize
from IPython.display import display
rng=np.random.default_rng(42)
plt.rcParams["figure.figsize"]=(7.2,4.5); plt.rcParams["axes.unicode_minus"]=False
pd.DataFrame({"package":["Python","NumPy","pandas","Matplotlib","SciPy"],"version":[platform.python_version(),np.__version__,pd.__version__,matplotlib.__version__,scipy.__version__]})

	package	version
0	Python	3.11.9
1	NumPy	1.26.4
2	pandas	2.2.2
3	Matplotlib	3.9.2
4	SciPy	1.13.1

架空データの作成

基準条件からの変化を $x_1=(T-190)/10, x_2=(P-80)/10$ と標準化します。品質応答には曲率・相互作用・測定ばらつきを含めます。

N=50
T=rng.uniform(175,210,N); P=rng.uniform(65,100,N); x1=(T-190)/10; x2=(P-80)/10
quality=96-2.8*(x1-.45)**2-1.9*(x2-.55)**2+.8*x1*x2-.2*x1**4+rng.normal(0,.45,N)
energy=18+2.4*x1+1.6*x2+.5*x1**2+rng.normal(0,.2,N)
trials=pd.DataFrame({"temperature_C":T,"pressure_MPa":P,"quality_score":quality,"energy_kWh":energy}).round(3)
display(trials.head()); display(trials.describe().round(2))
fig,ax=plt.subplots(1,2,figsize=(11,4.2)); s=ax[0].scatter(T,P,c=quality,cmap="viridis"); fig.colorbar(s,ax=ax[0],label="Quality score")
ax[0].set(title="Trial conditions and quality",xlabel="Temperature (°C)",ylabel="Pressure (MPa)"); ax[0].grid(True,alpha=.3)
ax[1].scatter(energy,quality); ax[1].set(title="Quality–energy trade-off",xlabel="Energy (kWh/lot)",ylabel="Quality score"); ax[1].grid(True,alpha=.3); plt.tight_layout(); plt.show()

	temperature_C	pressure_MPa	quality_score	energy_kWh
0	202.088	71.997	89.902	20.565
1	190.361	65.258	87.285	15.760
2	205.051	92.542	92.254	24.720
3	199.408	88.270	96.230	21.816
4	178.296	89.681	86.874	17.090

	temperature_C	pressure_MPa	quality_score	energy_kWh
count	50.00	50.00	50.00	50.00
mean	193.73	80.34	90.89	19.50
std	9.82	9.11	3.63	3.25
min	176.53	65.26	84.11	14.44
25%	186.05	72.13	87.76	17.03
50%	193.03	80.84	90.82	19.57
75%	202.20	88.23	94.11	21.77
max	209.15	98.67	96.23	26.60

No.041：非線形最適化とは何か

実務での意味

工程条件には適温があり、過不足の両方で品質が悪化します。曲がった応答面から品質とコストを両立する条件を選びます。

分析・モデル化の考え方

運転損失を品質損失とエネルギー負担の和、$f(x)=2.8(x_1-.45)^2+1.9(x_2-.55)^2-.8x_1x_2+.2x_1^4+.12(x_1+x_2)$ とします。

Pythonで確認する

def loss(x):
 a,b=x; return 2.8*(a-.45)**2+1.9*(b-.55)**2-.8*a*b+.2*a**4+.12*(a+b)
g=np.linspace(-1.5,2,150); X,Y=np.meshgrid(g,g); Z=np.vectorize(lambda a,b:loss([a,b]))(X,Y); r41=minimize(loss,[0,0],method="BFGS")
plt.contourf(190+10*X,80+10*Y,Z,25,cmap="viridis"); plt.colorbar(label="Operating loss"); plt.scatter(190+10*r41.x[0],80+10*r41.x[1],c="red",marker="*",s=180,label="Optimum")
plt.title("Nonlinear operating-loss surface"); plt.xlabel("Temperature (°C)"); plt.ylabel("Pressure (MPa)"); plt.grid(True,alpha=.3); plt.legend(); plt.tight_layout(); plt.show()
pd.Series({"temperature_C":190+10*r41.x[0],"pressure_MPa":80+10*r41.x[1],"minimum_loss":r41.fun,"success":r41.success})

temperature_C    194.998253
pressure_MPa      86.236484
minimum_loss      -0.084816
success                True
dtype: object

結果の読み取り

最適点は応答面の谷にあります。係数は不良費や電力費など同じ尺度へ換算し、根拠を明示します。

No.042：勾配とは何か

実務での意味

勾配は条件を少し動かしたときの損失感度で、負の勾配が局所的な改善方向です。

分析・モデル化の考え方

$\nabla f=(\partial f/\partial x_1,\partial f/\partial x_2)$。絶対値が大きい変数ほど調整効果が大きいと読めます。

Pythonで確認する

def grad(x):
 a,b=x; return np.array([5.6*(a-.45)-.8*b+.8*a**3+.12,3.8*(b-.55)-.8*a+.12])
pts=np.array([[-1,-.5],[0,0],[.8,.8],[1.5,1.2]])
pd.DataFrame([[190+10*a,80+10*b,loss([a,b]),*grad([a,b]),np.linalg.norm(grad([a,b]))] for a,b in pts],columns=["temperature_C","pressure_MPa","loss","grad_x1","grad_x2","gradient_norm"]).round(3)

	temperature_C	pressure_MPa	loss	grad_x1	grad_x2	gradient_norm
0	180.0	75.0	7.602	-8.40	-3.07	8.943
1	190.0	80.0	1.142	-2.40	-1.97	3.105
2	198.0	88.0	0.224	1.85	0.43	1.899
3	205.0	92.0	3.786	7.74	1.39	7.864

結果の読み取り

符号と大きさは地点で変わります。勾配ゼロは候補点であり、最小点かは曲率でも確認します。

No.043：勾配降下法

実務での意味

改善方向へ反復移動します。学習率が小さいと遅く、大きいと谷を飛び越えます。

分析・モデル化の考え方

更新式は $x^{(k+1)}=x^{(k)}-\alpha\nabla f(x^{(k)})$ です。停止には勾配ノルムを使います。

Pythonで確認する

def gd(x0,a,n=80):
 x=np.array(x0,dtype=float); h=[x.copy()]
 for _ in range(n):
  if np.linalg.norm(grad(x))<1e-7: break
  x-=a*grad(x); h.append(x.copy())
 return x,np.array(h)
rows=[]
for a in [.03,.12,.30]:
 x,h=gd([-1.2,1.5],a); v=[loss(z) for z in h]; plt.plot(v,label=f"alpha={a}"); rows.append([a,len(h)-1,v[-1],np.linalg.norm(grad(x))])
plt.title("Gradient descent convergence"); plt.xlabel("Iteration"); plt.ylabel("Operating loss"); plt.grid(True,alpha=.3); plt.legend(); plt.tight_layout(); plt.show()
pd.DataFrame(rows,columns=["learning_rate","iterations","final_loss","gradient_norm"]).round(6)

/var/folders/3y/fmw40k0x78xblvb3gkcyvy1h0000gn/T/ipykernel_25563/1166959326.py:2: RuntimeWarning: overflow encountered in scalar power
  a,b=x; return np.array([5.6*(a-.45)-.8*b+.8*a**3+.12,3.8*(b-.55)-.8*a+.12])
/var/folders/3y/fmw40k0x78xblvb3gkcyvy1h0000gn/T/ipykernel_25563/1702839514.py:5: RuntimeWarning: invalid value encountered in subtract
  x-=a*grad(x); h.append(x.copy())

	learning_rate	iterations	final_loss	gradient_norm
0	0.03	80	-0.084816	0.000134
1	0.12	30	-0.084816	0.000000
2	0.30	80	NaN	NaN

結果の読み取り

歩幅で収束速度が変わります。線探索、スケーリング、収束履歴保存が実務上重要です。

No.044：ニュートン法

実務での意味

傾きに加え曲率を使うため最適点近傍で高速ですが、2階微分の計算負荷に注意します。

分析・モデル化の考え方

Hessian $H$ を使い $x^{(k+1)}=x^{(k)}-H^{-1}\nabla f$。逆行列ではなく連立方程式を解きます。

Pythonで確認する

def hess(x): return np.array([[5.6+2.4*x[0]**2,-.8],[-.8,3.8]])
x=np.array([-1.2,1.5]); rows=[]
for k in range(15):
 rows.append([k,*x,loss(x),np.linalg.norm(grad(x))])
 if np.linalg.norm(grad(x))<1e-10: break
 x-=np.linalg.solve(hess(x),grad(x))
pd.DataFrame(rows,columns=["iteration","x1","x2","loss","gradient_norm"]).round(8)

	iteration	x1	x2	loss	gradient_norm
0	0	-1.200000	1.500000	11.228470	1.260723e+01
1	1	0.005619	0.519604	0.615376	2.784218e+00
2	2	0.518210	0.627518	-0.083794	1.112897e-01
3	3	0.499894	0.623662	-0.084816	4.123200e-04
4	4	0.499825	0.623647	-0.084816	1.000000e-08
5	5	0.499825	0.623647	-0.084816	0.000000e+00

結果の読み取り

少ない反復で勾配が急減します。比較は反復数でなく微分を含む総計算時間で行います。

No.045：準ニュートン法

実務での意味

Hessianを実装せず、勾配変化から曲率を近似して精度と負荷を両立します。

分析・モデル化の考え方

BFGSは移動量と勾配差からHessian近似を逐次更新します。

Pythonで確認する

rows=[]
for method in ["BFGS","CG"]:
 r=minimize(loss,[-1.2,1.5],jac=grad,method=method,options={"gtol":1e-9}); rows.append([method,r.success,r.nit,r.nfev,r.fun,190+10*r.x[0],80+10*r.x[1]])
pd.DataFrame(rows,columns=["method","success","iterations","function_evals","final_loss","temperature_C","pressure_MPa"]).round(6)

	method	success	iterations	function_evals	final_loss	temperature_C	pressure_MPa
0	BFGS	True	9	10	-0.084816	194.998255	86.236475
1	CG	True	8	19	-0.084816	194.998255	86.236475

結果の読み取り

成功フラグだけでなく、最終勾配と評価回数を記録します。

No.046：L-BFGS

実務での意味

多数の変数を持つ問題では、直近履歴だけを持つL-BFGSがメモリを抑えます。

分析・モデル化の考え方

100工程の個別目標との差と隣接工程の急変を同時に罰します。L-BFGS-Bは上下限も扱えます。

Pythonで確認する

q=100; target=.5+.25*np.sin(np.linspace(0,4*np.pi,q))
def lf(z): return np.sum((z-target)**2)+.4*np.sum(np.diff(z)**2)+.02*np.sum(z**4)
def lg(z):
 g=2*(z-target)+.08*z**3; d=np.diff(z); g[:-1]-=.8*d; g[1:]+=.8*d; return g
r46=minimize(lf,np.zeros(q),jac=lg,method="L-BFGS-B",bounds=[(-1,1)]*q)
plt.plot(target,"--",label="Individual target"); plt.plot(r46.x,label="Optimized setting"); plt.title("L-BFGS-B for 100 process settings"); plt.xlabel("Process index"); plt.ylabel("Scaled setting"); plt.grid(True,alpha=.3); plt.legend(); plt.tight_layout(); plt.show()
pd.Series({"success":r46.success,"iterations":r46.nit,"final_loss":r46.fun})

success           True
iterations           8
final_loss    0.231526
dtype: object

結果の読み取り

個別目標を追いつつ急変を抑えました。メモリ、計算時間、再計算頻度も要件化します。

No.047：ラグランジュ未定乗数法

実務での意味

設備負荷合計が一定なら、制約境界上の最良点を求めます。

分析・モデル化の考え方

$q=(x_1-1)^2+2(x_2-1)^2$, $h=x_1+x_2-1=0$ として $\mathcal L=q+\lambda h$ の停留条件を解きます。

Pythonで確認する

A=np.array([[2,0,1],[0,4,1],[1,1,0]],float); b=np.array([2,4,1],float); a,bv,lam=np.linalg.solve(A,b)
pd.Series({"x1":a,"x2":bv,"lambda":lam,"constraint_residual":a+bv-1,"quality_loss":(a-1)**2+2*(bv-1)**2}).round(6)

x1                     0.333333
x2                     0.666667
lambda                 1.333333
constraint_residual    0.000000
quality_loss           0.666667
dtype: float64

結果の読み取り

乗数は制約緩和の局所価値です。符号は定義に依存するため、右辺変更時の目的値変化で確認します。

No.048：KKT条件

実務での意味

どの制約が解を縛るかを特定し、設備増強や規格見直しの材料にします。

分析・モデル化の考え方

$g_i(x)\le0$ の下で、主・双対実行可能性、停留条件、相補性 $\mu_i g_i=0$ を点検します。

Pythonで確認する

def ko(x): return (x[0]-1)**2+(x[1]-1)**2
r48=minimize(ko,[.2,.2],method="SLSQP",constraints=[{"type":"ineq","fun":lambda x:1.2-x.sum()}]); x=r48.x; mu=2*(1-x[0]); gv=x.sum()-1.2
st=np.array([2*(x[0]-1)+mu,2*(x[1]-1)+mu])
pd.Series({"x1":x[0],"x2":x[1],"g(x)<=0":gv,"mu>=0":mu,"stationarity_norm":np.linalg.norm(st),"complementarity":mu*gv}).round(8)

x1                   0.6
x2                   0.6
g(x)<=0              0.0
mu>=0                0.8
stationarity_norm    0.0
complementarity      0.0
dtype: float64

結果の読み取り

制約は活性でKKT残差もほぼゼロです。活性制約と乗数は追加能力の価値を示します。

No.049：内点法

実務での意味

不等式制約の内側を進み、境界上の解へ近づきます。

分析・モデル化の考え方

$x_1,x_2>0$, $x_1+x_2<1.2$ に対して対数バリアを加え、徐々に弱めます。

Pythonで確認する

def barrier(x,t):
 s=1.2-x.sum()
 return np.inf if min(x)<=0 or s<=0 else ko(x)-(np.log(x[0])+np.log(x[1])+np.log(s))/t
x=np.array([.3,.3]); path=[]
for t in [1,3,10,30,100,300,1000]:
 r=minimize(lambda z:barrier(z,t),x,method="Nelder-Mead"); x=r.x; path.append([t,*x,1.2-x.sum(),ko(x)])
df=pd.DataFrame(path,columns=["t","x1","x2","slack","objective"]); display(df.round(6))
plt.plot(df.x1,df.x2,marker="o"); plt.scatter(.6,.6,c="red",marker="*",s=160,label="Optimum"); plt.title("Central path of barrier method"); plt.xlabel("x1"); plt.ylabel("x2"); plt.grid(True,alpha=.3); plt.legend(); plt.tight_layout(); plt.show()

	t	x1	x2	slack	objective
0	1	0.449621	0.449651	0.300727	0.605800
1	3	0.499973	0.500000	0.200026	0.500026
2	10	0.553489	0.553382	0.093129	0.398840
3	30	0.581397	0.581342	0.037261	0.350503
4	100	0.594008	0.593921	0.012071	0.329730
5	300	0.597956	0.597923	0.004121	0.323306
6	1000	0.599346	0.599408	0.001245	0.320997

結果の読み取り

バリアを弱めると制約余裕が減ります。許容誤差を物理単位へ戻し安全余裕を確認します。

No.050：非線形最適化と機械学習

実務での意味

予測モデルを次に運転する条件という意思決定へつなぎます。

分析・モデル化の考え方

試作から非線形回帰 $\hat y(x)$ を学習し、$\min\{-\hat y(x)+\rho E(x)\}$ を試作範囲内で解きます。

Pythonで確認する

def feat(a,b):
 a,b=np.asarray(a),np.asarray(b); return np.column_stack([np.ones_like(a),a,b,a*a,a*b,b*b,a**3,a**4])
coef,*_=np.linalg.lstsq(feat(x1,x2),quality,rcond=None); pred=feat(x1,x2)@coef; rmse=np.sqrt(np.mean((quality-pred)**2))
def pq(x): return float((feat([x[0]],[x[1]])@coef).item())
def dl(x): return -pq(x)+.12*(18+2.4*x[0]+1.6*x[1]+.5*x[0]**2)
r50=minimize(dl,[0,0],method="L-BFGS-B",bounds=[(x1.min(),x1.max()),(x2.min(),x2.max())])
display(pd.Series({"training_RMSE":rmse,"recommended_temperature_C":190+10*r50.x[0],"recommended_pressure_MPa":80+10*r50.x[1],"predicted_quality":pq(r50.x),"success":r50.success}).round(4))
plt.scatter(pred,quality); lo=min(pred.min(),quality.min()); hi=max(pred.max(),quality.max()); plt.plot([lo,hi],[lo,hi],"--",color="black",label="Ideal"); plt.title("Observed vs predicted quality"); plt.xlabel("Predicted quality score"); plt.ylabel("Observed quality score"); plt.grid(True,alpha=.3); plt.legend(); plt.tight_layout(); plt.show()

training_RMSE                   0.43307
recommended_temperature_C    194.965778
recommended_pressure_MPa      86.175973
predicted_quality             96.258864
success                            True
dtype: object

結果の読み取り

学習内RMSEだけでは不十分です。検証データ、予測区間、外挿防止、推奨条件での確認試験が必要です。

対象ノックを通して見える実務上の示唆

• 目的係数には不良費・電力費・納期影響などの根拠が必要です。
• 変数のスケーリングは収束と許容誤差の管理に効きます。
• KKT乗数と活性制約は設備増強の優先順位に使えます。
• 初期値、収束履歴、制約残差、再計算条件を保存します。
• 予測精度と最適化結果は分けて検証し、外挿を防ぎます。

実務導入する場合に必要なこと

論点	確認事項	成果物
KPI	品質・原価・納期の換算	目的関数定義書
制約	設備限界、規格、安全余裕	制約一覧・根拠
データ	測定系、欠測、試験範囲	データ辞書
検証	現行比較、感度、外挿	確認試験計画
運用	承認、再計算、異常時対応	SOP・監視指標

まず人が推奨値を確認する意思決定支援として始め、効果と安全性を検証しながら自動化範囲を広げます。

まとめ

非線形応答の定式化から、勾配・曲率を使う解法、制約付き解の確認、機械学習による条件推薦までを確認しました。最適化は計算だけでなく、目的・制約を明文化し改善余地とリスクを共有する意思決定プロセスです。

法人向けのご相談

数理工房では、製造条件最適化、生産計画、シミュレーション、予測と最適化の統合について、課題整理・PoC・運用設計をご支援します。

📩 お問い合わせ: surikobo.co.jp/contact
まずはお気軽にご相談ください。