品質と設備制約を同時に守る：凸最適化で成形条件を説明可能に決める10本ノック

架空の樹脂成形工程を題材に、品質損失・電力負荷・設備能力を同じモデルで扱い、「なぜその条件なのか」を現場へ説明できる条件設計を学びます。対象は No.051〜No.060（凸最適化） です。

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

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

成形品質は、樹脂温度と保持圧力を上げれば単調に改善するとは限りません。条件不足では充填不良、条件過多では熱劣化やバリが増え、さらに設備能力、電力、標準条件からの変更幅も守る必要があります。本記事では、これらを「候補条件を勘で試す問題」から「目的と制約を明記して再計算できる問題」へ変換します。

現場でよくある状況

• 熟練者ごとに推奨条件が異なり、根拠が口頭でしか残っていない
• 品質、サイクルタイム、電力のKPIが別々に管理され、全体最適にならない
• 設備上限ぎりぎりの条件が出ても、上限緩和に投資する価値を説明できない
• センサーを増やした結果、説明変数が増えすぎてモデルの運用が不安定になる

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

目的が複数あり、変数同士が連動し、制約も同時に満たす必要があるためです。ただし、実行可能領域が凸集合、目的関数が凸関数なら、局所解が大域解となり、最適性条件や双対変数を使って結果を監査できます。一方、実工程に段取りの有無や複数の安定領域が含まれると凸性は崩れます。本記事では「凸最適化が使える範囲」と「使えない兆候」の両方を扱います。

今回扱うノックの全体像

No.	テーマ	製造業での判断
051	凸集合	条件を平均しても安全か
052	凸関数	中間条件の損失を上から評価できるか
053	凸最適化の重要性	初期値に依存せず再現できるか
054	一次条件	解が最適であることを検算できるか
055	ラグランジュ双対	制約から下界と経済価値を得られるか
056	強双対性	主問題と双対問題が一致するか
057	制約付き凸最適化	実設備の範囲内で条件を決められるか
058	正則化	変更量・複雑さを抑えられるか
059	Lasso	重要なセンサーを絞れるか
060	凸最適化の限界	非凸性を見落としていないか

Python 環境の準備

外部データは使いません。乱数生成器を固定し、同じデータと結果を再現できるようにします。最適化にはSciPy、Lassoにはscikit-learn、可視化にはMatplotlibを使います。

import platform
import numpy as np
import pandas as pd
import matplotlib
import matplotlib.pyplot as plt
import scipy
import sklearn
import japanize_matplotlib

from IPython.display import display
from scipy.optimize import minimize
from sklearn.linear_model import Lasso
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import StandardScaler

SEED = 20260712
rng = np.random.default_rng(SEED)
plt.rcParams["figure.figsize"] = (7.2, 4.5)
plt.rcParams["axes.unicode_minus"] = False

pd.DataFrame({
    "package": ["Python", "NumPy", "pandas", "SciPy", "scikit-learn", "Matplotlib"],
    "version": [platform.python_version(), np.__version__, pd.__version__, scipy.__version__, sklearn.__version__, matplotlib.__version__],
})

	package	version
0	Python	3.13.1
1	NumPy	2.5.1
2	pandas	3.0.3
3	SciPy	1.18.0
4	scikit-learn	1.9.0
5	Matplotlib	3.11.0

架空データの作成

対象は樹脂成形工程です。温度 $T$ と保持圧力 $P$ を扱いやすい無次元量

$x_1=\frac{T-190}{10},\qquad x_2=\frac{P-80}{10}$

に変換します。標準化後の品質損失を

$f(x)=2(x_1-0.6)^2+1.5(x_2-0.4)^2+0.6(x_1-x_2)^2$

と置きます。二次形式のHessianは正定値なので、この損失は狭義凸です。観測データには測定ばらつきを加えますが、最適化に使う基準モデル自体は明示的に管理します。

n = 90
x1_obs = rng.uniform(-1.2, 1.5, n)
x2_obs = rng.uniform(-1.0, 1.4, n)

def quality_loss(x):
    x = np.asarray(x)
    return 2.0 * (x[..., 0] - 0.6) ** 2 + 1.5 * (x[..., 1] - 0.4) ** 2 + 0.6 * (x[..., 0] - x[..., 1]) ** 2

X_obs = np.column_stack([x1_obs, x2_obs])
loss_obs = quality_loss(X_obs) + rng.normal(0, 0.18, n)
process_df = pd.DataFrame({
    "樹脂温度_C": 190 + 10 * x1_obs,
    "保持圧力_MPa": 80 + 10 * x2_obs,
    "品質損失_指数": loss_obs,
})
display(process_df.head())
display(process_df.describe().round(2))

	樹脂温度_C	保持圧力_MPa	品質損失_指数
0	196.697726	81.129631	-0.211684
1	191.024661	77.871133	0.990317
2	201.472181	90.085739	0.922645
3	189.287092	92.945395	3.540957
4	178.914627	91.178000	9.526144

	樹脂温度_C	保持圧力_MPa	品質損失_指数
count	90.00	90.00	90.00
mean	191.54	82.06	3.06
std	8.15	6.57	2.43
min	178.16	70.52	-0.21
25%	184.73	76.48	1.31
50%	191.38	82.21	2.18
75%	198.36	87.04	4.83
max	204.94	92.97	9.53

No.051：凸集合とは何か

実務での意味

二つの安全な運転条件を結ぶ途中の条件もすべて安全なら、その条件領域は凸集合です。現場で条件を徐々に切り替える場合や、複数製品の平均負荷を考える場合に重要です。

分析・モデル化の考え方

集合 $C$ が凸であるとは、任意の $x,y\in C$ と $0\le\theta\le1$ に対し、$\theta x+(1-\theta)y\in C$ が成立することです。ここでは温度・圧力の上下限と合計負荷 $x_1+x_2\le1.2$ の共通部分を実行可能領域とします。線形不等式の共通部分は凸集合です。

Pythonで確認する

def feasible(x):
    a, b = np.asarray(x)
    return (-1 <= a <= 1.4) and (-1 <= b <= 1.3) and (a + b <= 1.2)

A, B = np.array([-0.8, 1.1]), np.array([1.1, -0.4])
theta = np.linspace(0, 1, 21)
segment = theta[:, None] * A + (1 - theta[:, None]) * B

g = np.linspace(-1.3, 1.6, 240)
G1, G2 = np.meshgrid(g, g)
mask = (G1 >= -1) & (G1 <= 1.4) & (G2 >= -1) & (G2 <= 1.3) & (G1 + G2 <= 1.2)
plt.contourf(190 + 10 * G1, 80 + 10 * G2, mask.astype(int), levels=[0.5, 1.5], alpha=0.35, colors=["tab:blue"])
plt.plot(190 + 10 * segment[:, 0], 80 + 10 * segment[:, 1], "o-", color="tab:orange", label="2条件を結ぶ線分")
plt.scatter(190 + 10 * np.array([A[0], B[0]]), 80 + 10 * np.array([A[1], B[1]]), color="black", zorder=3, label="端点")
plt.title("凸な実行可能領域と条件間の補間")
plt.xlabel("樹脂温度 [°C]"); plt.ylabel("保持圧力 [MPa]")
plt.grid(True, alpha=0.3); plt.legend(); plt.tight_layout(); plt.show()
print(f"線分上の全点が実行可能: {all(map(feasible, segment))}")

線分上の全点が実行可能: True

結果の読み取り

青い領域内の二点を結ぶ橙色の線分は、すべて領域内に残ります。この性質があれば、条件補間によって突然制約違反になる心配を減らせます。ただし、温度帯の中央だけ材料が不安定になるような「穴のある領域」は非凸です。安全領域の定義は設備仕様だけでなく、品質保証部門の知見でも検証します。

No.052：凸関数とは何か

実務での意味

凸な損失関数では、極端な二条件の中間を採用したときの損失が、両端損失の加重平均を超えません。条件変更のリスクを評価しやすく、改善方向も一貫します。

分析・モデル化の考え方

$f$ が凸であるとは、$f(\theta x+(1-\theta)y)\le\theta f(x)+(1-\theta)f(y)$ が成立することです。滑らかな二次関数ではHessianの固有値がすべて非負なら凸です。本モデルの曲率を数値で確認し、Jensenの不等式の余裕も可視化します。

Pythonで確認する

H = np.array([[5.2, -1.2], [-1.2, 4.2]])
eig = np.linalg.eigvalsh(H)
lhs = quality_loss(segment)
rhs = theta * quality_loss(A) + (1 - theta) * quality_loss(B)
display(pd.DataFrame({"Hessian固有値": eig}).round(3))
plt.plot(theta, lhs, marker="o", label="中間条件の実損失")
plt.plot(theta, rhs, "--", label="端点損失の加重平均（上界）")
plt.title("凸関数が満たすJensenの不等式")
plt.xlabel("端点Aの混合比 θ"); plt.ylabel("品質損失指数")
plt.grid(True, alpha=0.3); plt.legend(); plt.tight_layout(); plt.show()
print(f"Jensen不等式の最大違反量: {np.max(lhs-rhs):.3e}")

	Hessian固有値
0	3.4
1	6.0

Jensen不等式の最大違反量: 0.000e+00

結果の読み取り

Hessianの固有値がともに正で、実損失は上界以下です。したがって本モデルは狭義凸で、谷底は一つです。実務では、推定した応答曲面が凸に見えても、観測範囲外まで保証されたわけではありません。適用範囲と曲率の不確実性をモデル台帳に残します。

No.053：凸最適化がなぜ重要なのか

実務での意味

担当者や初期条件が変わっても同じ推奨条件へ到達できれば、条件表の再現性と監査性が高まります。凸問題では局所最適解が大域最適解なので、探索開始点への依存を大幅に減らせます。

分析・モデル化の考え方

異なる初期値から同じ凸目的関数を最小化し、解と目的値のばらつきを確認します。計算上の許容誤差は残るため、完全一致ではなく物理単位で十分小さい差かを評価します。

Pythonで確認する

starts = np.array([[-1, -1], [-1, 1], [1.3, -0.8], [1.2, 1.1], [0, 0]])
rows = []
for s in starts:
    result = minimize(quality_loss, s, method="BFGS", options={"gtol": 1e-10})
    rows.append([*s, *result.x, result.fun, result.nit, result.success])
convex_runs = pd.DataFrame(rows, columns=["初期x1", "初期x2", "解x1", "解x2", "最終損失", "反復数", "成功"])
convex_runs["最適温度_C"] = 190 + 10 * convex_runs["解x1"]
convex_runs["最適圧力_MPa"] = 80 + 10 * convex_runs["解x2"]
display(convex_runs.round(6))
print(f"初期値間の最大解差: {np.ptp(convex_runs[['解x1','解x2']].to_numpy(), axis=0).max():.2e}")

	初期x1	初期x2	解x1	解x2	最終損失	反復数	成功	最適温度_C	最適圧力_MPa
0	-1.0	-1.0	0.564706	0.447059	0.014118	5	False	195.647059	84.470588
1	-1.0	1.0	0.564706	0.447059	0.014118	6	True	195.647059	84.470588
2	1.3	-0.8	0.564706	0.447059	0.014118	4	False	195.647059	84.470588
3	1.2	1.1	0.564706	0.447059	0.014118	4	True	195.647059	84.470588
4	0.0	0.0	0.564706	0.447059	0.014118	6	True	195.647059	84.470588

初期値間の最大解差: 1.45e-10

結果の読み取り

すべての開始点がほぼ同じ条件と損失へ収束します。これは「ソルバーが成功と表示した」以上の再現性確認です。本番では、目的値、制約残差、停止理由、入力データ版、ソルバー版も保存し、再計算可能にします。

No.054：一次条件と最適性条件

実務での意味

推奨条件が単なる探索結果ではなく、微小な条件変更では改善できない点だと検算します。条件変更会議で「少し温度を上げれば良いのでは」という提案に定量的に答えられます。

分析・モデル化の考え方

微分可能な凸関数では、制約なしの点 $x^*$ が最適である必要十分条件は $\nabla f(x^*)=0$ です。また一次の支持超平面

$f(y)\ge f(x)+\nabla f(x)^\top(y-x)$

が任意の $x,y$ で成立します。勾配ノルムと一次下界を確認します。

Pythonで確認する

def quality_grad(x):
    a, b = x
    return np.array([4 * (a - 0.6) + 1.2 * (a - b), 3 * (b - 0.4) - 1.2 * (a - b)])

x_star = minimize(quality_loss, [0, 0], jac=quality_grad, method="BFGS", tol=1e-12).x
x_base = np.array([-0.3, 0.8])
y_test = rng.uniform(-1, 1.3, size=(200, 2))
first_order_lower = quality_loss(x_base) + (y_test - x_base) @ quality_grad(x_base)
gap = quality_loss(y_test) - first_order_lower
display(pd.Series({
    "最適x1": x_star[0], "最適x2": x_star[1],
    "勾配ノルム": np.linalg.norm(quality_grad(x_star)),
    "一次下界の最小余裕": gap.min(),
}).round(10))

最適x1         0.564706
最適x2         0.447059
勾配ノルム        0.000000
一次下界の最小余裕    0.024987
dtype: float64

結果の読み取り

最適点の勾配ノルムはほぼゼロで、試した点すべてで一次下界が実際の損失以下です。制約がある場合は勾配ゼロにならないことがあり、そのときは制約の法線と目的勾配の釣り合い、すなわちKKT条件で検算します。

No.055：双対性とラグランジュ双対

実務での意味

設備負荷の上限が最適条件を縛る場合、上限を少し緩める価値を金額や損失指数で表せます。これが増強投資やボトルネック改善の優先順位につながります。

分析・モデル化の考え方

説明を明確にするため、目標条件 $c=(0.9,0.7)$ への二乗距離を最小化し、能力制約 $x_1+x_2\le b$ を課します。

$\min_x\ \frac12\lVert x-c\rVert^2,\qquad x_1+x_2-b\le0$

ラグランジアン $L(x,\lambda)=\frac12\lVert x-c\rVert^2+\lambda(x_1+x_2-b)$ を $x$ について最小化すると、双対関数は $g(\lambda)=-\lambda^2+\lambda(c_1+c_2-b)$ です。任意の $\lambda\ge0$ で $g(\lambda)$ は主問題最適値の下界になります。

Pythonで確認する

c = np.array([0.9, 0.7]); capacity = 1.2
dual_lambda = np.linspace(0, 0.8, 161)
dual_value = -dual_lambda**2 + dual_lambda * (c.sum() - capacity)
lambda_star = (c.sum() - capacity) / 2
primal_x = c - lambda_star * np.ones(2)
primal_value = 0.5 * np.sum((primal_x - c)**2)
plt.plot(dual_lambda, dual_value, label="双対関数 g(λ)")
plt.axhline(primal_value, color="tab:red", linestyle="--", label="主問題の最適値")
plt.scatter([lambda_star], [dual_value.max()], color="black", zorder=3)
plt.title("双対関数による最適値の下界")
plt.xlabel("双対変数 λ"); plt.ylabel("目的値／下界")
plt.grid(True, alpha=0.3); plt.legend(); plt.tight_layout(); plt.show()
display(pd.Series({"最適λ": lambda_star, "最適x1": primal_x[0], "最適x2": primal_x[1], "主問題値": primal_value, "双対最大値": dual_value.max()}).round(4))

最適λ      0.20
最適x1     0.70
最適x2     0.50
主問題値     0.04
双対最大値    0.04
dtype: float64

結果の読み取り

双対関数はどの非負のλでも主問題値を超えず、最適λで一致します。λは能力上限をわずかに緩和したときの目的値改善率です。ただし、尺度は目的関数の単位に依存します。投資判断に使うなら品質損失指数を円へ換算し、有限差分でも符号と大きさを検証します。

No.056：強双対性

実務での意味

主問題の最良値と双対問題の下界が一致すれば、解の品質を数値的に証明できます。計算時間を打ち切る大規模問題でも、双対ギャップは「あとどれだけ改善余地があるか」の指標になります。

分析・モデル化の考え方

凸問題で、すべての不等式制約を厳密に満たす点が存在するというSlater条件などが成立すると強双対性が得られます。前問では $(0,0)$ が $x_1+x_2<1.2$ を満たすため、この条件が成立します。能力上限を変えて主値、双対値、感度を比較します。

Pythonで確認する

records = []
for b in [0.9, 1.0, 1.1, 1.2, 1.3, 1.4, 1.6, 1.8]:
    lam = max(0.0, (c.sum() - b) / 2)
    x = c - lam * np.ones(2)
    pv = 0.5 * np.sum((x - c)**2)
    dv = -lam**2 + lam * (c.sum() - b)
    records.append([b, *x, lam, pv, dv, pv-dv])
duality_df = pd.DataFrame(records, columns=["能力上限b", "x1", "x2", "λ", "主問題値", "双対値", "双対ギャップ"])
display(duality_df.round(8))
plt.plot(duality_df["能力上限b"], duality_df["主問題値"], "o-", label="最小損失")
plt.plot(duality_df["能力上限b"], duality_df["λ"], "s--", label="制約緩和の限界価値 λ")
plt.title("設備能力と損失・双対変数の関係")
plt.xlabel("能力上限 b"); plt.ylabel("目的値／双対変数")
plt.grid(True, alpha=0.3); plt.legend(); plt.tight_layout(); plt.show()

	能力上限b	x1	x2	λ	主問題値	双対値	双対ギャップ
0	0.9	0.55	0.35	0.35	0.1225	0.1225	-0.0
1	1.0	0.60	0.40	0.30	0.0900	0.0900	0.0
2	1.1	0.65	0.45	0.25	0.0625	0.0625	0.0
3	1.2	0.70	0.50	0.20	0.0400	0.0400	0.0
4	1.3	0.75	0.55	0.15	0.0225	0.0225	0.0
5	1.4	0.80	0.60	0.10	0.0100	0.0100	0.0
6	1.6	0.90	0.70	0.00	0.0000	0.0000	0.0
7	1.8	0.90	0.70	0.00	0.0000	-0.0000	0.0

結果の読み取り

すべての能力水準で双対ギャップは丸め誤差の範囲です。能力が厳しいほどλが大きく、増強の限界価値が高くなります。上限が目標条件の合計1.6以上になると制約は非活性となり、λはゼロです。実務では一度の増強幅が大きい場合、局所感度λだけでなく上限変更後に再最適化します。

No.057：制約付き凸最適化

実務での意味

品質だけを追った理想条件が設備範囲外なら実行できません。温度・圧力の上下限、合計負荷、標準条件からの変更幅を同時に満たす実装可能な条件を求めます。

分析・モデル化の考え方

凸目的関数に線形不等式制約を組み合わせます。SciPyのSLSQPで解きますが、結果は成功フラグだけでなく、全制約の余裕と、近傍候補との目的値比較でも点検します。制約は $g(x)\ge0$ の形で実装します。

Pythonで確認する

constraints = [
    {"type": "ineq", "fun": lambda x: x[0] + 1.0},
    {"type": "ineq", "fun": lambda x: 1.4 - x[0]},
    {"type": "ineq", "fun": lambda x: x[1] + 1.0},
    {"type": "ineq", "fun": lambda x: 1.3 - x[1]},
    {"type": "ineq", "fun": lambda x: 0.75 - x[0] - x[1]},
]
r57 = minimize(quality_loss, [0, 0], jac=quality_grad, method="SLSQP", constraints=constraints,
               options={"ftol": 1e-12, "maxiter": 500})
margins = np.array([con["fun"](r57.x) for con in constraints])
display(pd.Series({
    "solver_success": r57.success, "最適温度_C": 190 + 10*r57.x[0],
    "最適圧力_MPa": 80 + 10*r57.x[1], "品質損失": r57.fun,
    "最小制約余裕": margins.min(), "反復数": r57.nit,
}))

Z = quality_loss(np.stack([G1, G2], axis=-1))
plt.contour(190+10*G1, 80+10*G2, Z, levels=12, cmap="viridis")
mask57 = (G1 >= -1) & (G1 <= 1.4) & (G2 >= -1) & (G2 <= 1.3) & (G1 + G2 <= 0.75)
plt.contourf(190+10*G1, 80+10*G2, mask57.astype(int), levels=[0.5, 1.5], alpha=0.18, colors=["tab:blue"])
plt.scatter(190+10*r57.x[0], 80+10*r57.x[1], marker="*", s=180, color="tab:red", label="制約付き最適解")
plt.title("実行可能領域内の品質損失最小化")
plt.xlabel("樹脂温度 [°C]"); plt.ylabel("保持圧力 [MPa]")
plt.grid(True, alpha=0.3); plt.legend(); plt.tight_layout(); plt.show()

solver_success          True
最適温度_C            194.449153
最適圧力_MPa           83.050847
品質損失                0.073347
最小制約余裕                   0.0
反復数                        5
dtype: object

結果の読み取り

最適解は合計負荷の境界上にあり、設備能力が品質改善を制限しています。最小制約余裕が負でないため数値上は実行可能です。ただし、制御誤差を考えると境界値をそのまま設定値にしてはいけません。実運用では温度・圧力のばらつきに応じた安全余裕を追加し、その条件で再最適化します。

No.058：正則化と凸最適化

実務での意味

毎回大きく条件を変える提案は、品質が良くても段取り負荷や現場の受容性を悪化させます。正則化は、性能と「変更の少なさ」「運用の単純さ」を同じ目的関数で調整します。

分析・モデル化の考え方

前回設定 $x_{prev}=(-0.2,0.1)$ からの変更にL2罰則を加えます。

$\min_x\ f(x)+\alpha\lVert x-x_{prev}\rVert_2^2$

$\alpha$ が大きいほど現状維持を優先します。いずれも凸関数なので、和も凸です。複数の $\alpha$ で品質損失と変更量のトレードオフを比較します。

Pythonで確認する

x_prev = np.array([-0.2, 0.1])
reg_rows = []
for alpha in [0, 0.1, 0.3, 1, 3, 10]:
    objective = lambda x, a=alpha: quality_loss(x) + a*np.sum((x-x_prev)**2)
    r = minimize(objective, x_prev, method="BFGS")
    reg_rows.append([alpha, *r.x, quality_loss(r.x), np.linalg.norm(r.x-x_prev), r.fun])
reg_df = pd.DataFrame(reg_rows, columns=["正則化強度α", "x1", "x2", "品質損失", "設定変更量", "正則化込み目的値"])
display(reg_df.round(4))
plt.plot(reg_df["設定変更量"], reg_df["品質損失"], "o-")
for _, row in reg_df.iterrows():
    plt.annotate(f"α={row['正則化強度α']:g}", (row["設定変更量"], row["品質損失"]), xytext=(4, 4), textcoords="offset points")
plt.title("品質改善と設定変更量のトレードオフ")
plt.xlabel("前回設定からの変更量"); plt.ylabel("品質損失指数")
plt.grid(True, alpha=0.3); plt.tight_layout(); plt.show()

	正則化強度α	x1	x2	品質損失	設定変更量	正則化込み目的値
0	0.0	0.5647	0.4471	0.0141	0.8398	0.0141
1	0.1	0.5308	0.4220	0.0174	0.7986	0.0812
2	0.3	0.4718	0.3805	0.0384	0.7280	0.1974
3	1.0	0.3259	0.2889	0.1696	0.5588	0.4819
4	3.0	0.1277	0.1915	0.5139	0.3402	0.8611
5	10.0	-0.0573	0.1294	0.9949	0.1457	1.2071

結果の読み取り

正則化を強めると設定変更量は減りますが、品質損失は増えます。αは数学的に自動決定される経営判断ではありません。条件変更に必要な停止時間、検証費用、作業リスクをできれば円換算し、操業方針として合意します。定期改定と緊急改定でαを分ける設計も有効です。

No.059：Lassoと凸最適化

実務での意味

多数のセンサーから品質を説明すると、不要な信号までモデルへ入り、保守費用や誤警報が増えます。Lassoは予測誤差を抑えながら一部係数をゼロにし、監視対象を絞る候補を作ります。

分析・モデル化の考え方

Lassoは

$\min_{\beta_0,\beta}\ \frac{1}{2n}\lVert y-\beta_0-X\beta\rVert_2^2+\alpha\lVert\beta\rVert_1$

を解きます。二乗誤差とL1ノルムはいずれも凸です。単位差による不公平を避けるため説明変数を標準化します。ここでは真に影響する変数を温度、圧力、冷却時間、材料水分の4つに限定した架空データを生成します。

Pythonで確認する

n_lasso, p = 160, 12
sensor_names = ["樹脂温度", "保持圧力", "冷却時間", "材料水分", "金型温度", "射出速度", "背圧", "外気温", "油温", "振動", "待機時間", "ロット番号"]
X_sensor = rng.normal(size=(n_lasso, p))
true_beta = np.array([2.4, -1.8, 1.2, -0.9] + [0.0]*(p-4))
y_quality = 10 + X_sensor @ true_beta + rng.normal(0, 0.9, n_lasso)

lasso_rows = []
coef_paths = []
alphas = np.logspace(-2.2, 0.2, 20)
for alpha in alphas:
    model = make_pipeline(StandardScaler(), Lasso(alpha=alpha, max_iter=20000))
    model.fit(X_sensor, y_quality)
    coef = model[-1].coef_
    pred = model.predict(X_sensor)
    lasso_rows.append([alpha, np.count_nonzero(np.abs(coef) > 1e-8), np.sqrt(np.mean((y_quality-pred)**2))])
    coef_paths.append(coef)
lasso_df = pd.DataFrame(lasso_rows, columns=["α", "非ゼロ係数数", "学習RMSE"])
selected_idx = np.argmin(np.abs(alphas-0.18))
coef_table = pd.DataFrame({"センサー": sensor_names, "Lasso係数": coef_paths[selected_idx], "真の係数（生成時）": true_beta})
display(lasso_df.iloc[[0, 5, 10, 15, 19]].round(3))
display(coef_table.reindex(coef_table["Lasso係数"].abs().sort_values(ascending=False).index).round(3))

coef_paths = np.asarray(coef_paths)
for j, name in enumerate(sensor_names):
    plt.plot(alphas, coef_paths[:, j], label=name)
plt.xscale("log")
plt.title("Lassoの正則化パス")
plt.xlabel("正則化強度 α（対数軸）"); plt.ylabel("標準化後の回帰係数")
plt.grid(True, alpha=0.3); plt.legend(bbox_to_anchor=(1.02, 1), loc="upper left", fontsize=8); plt.tight_layout(); plt.show()

	α	非ゼロ係数数	学習RMSE
0	0.006	12	0.893
5	0.027	8	0.898
10	0.116	4	0.940
15	0.495	4	1.452
19	1.585	1	2.778

	センサー	Lasso係数	真の係数（生成時）
0	樹脂温度	2.074	2.4
1	保持圧力	-1.513	-1.8
2	冷却時間	1.043	1.2
3	材料水分	-0.766	-0.9
4	金型温度	-0.000	0.0
5	射出速度	0.000	0.0
6	背圧	0.000	0.0
7	外気温	0.000	0.0
8	油温	-0.000	0.0
9	振動	-0.000	0.0
10	待機時間	0.000	0.0
11	ロット番号	0.000	0.0

結果の読み取り

αを強めるほど非ゼロ係数が減り、主要4変数が相対的に残ります。ただし、Lassoの選択は因果関係の証明ではありません。相関の強いセンサー群では代表が入れ替わり、学習データだけのRMSEは楽観的です。時系列順の検証、工程知識、安全上必須の監視項目、センサー故障時の代替性を合わせて選びます。

No.060：凸最適化の限界

実務での意味

凸モデルは解きやすく説明もしやすい一方、工程に複数の安定運転帯、段取りの有無、整数台数、起動停止などがあると現実を表し切れません。「解けるモデル」に現場を合わせると、見かけ上きれいでも実行不能な提案になります。

分析・モデル化の考え方

非凸関数には複数の局所解があり、一次条件を満たしても大域最適とは限りません。例として、一変数の損失

$h(z)=(z^2-1)^2+0.18z$

を考えます。二階微分 $h''(z)=12z^2-4$ は原点付近で負なので非凸です。異なる初期値から局所最適化し、全域グリッドの最良値と比較します。

Pythonで確認する

def nonconvex_loss(z):
    z = np.asarray(z)
    return (z**2 - 1)**2 + 0.18*z

z_grid = np.linspace(-1.8, 1.8, 1200)
nc_rows = []
for start in [-1.6, -0.5, 0.2, 0.8, 1.6]:
    r = minimize(lambda q: float(nonconvex_loss(q[0])), [start], method="BFGS")
    nc_rows.append([start, r.x[0], r.fun, r.success])
nc_df = pd.DataFrame(nc_rows, columns=["初期値", "到達点", "局所目的値", "成功"])
global_i = np.argmin(nonconvex_loss(z_grid))
display(nc_df.round(5))
print(f"グリッドで確認した大域最良候補: z={z_grid[global_i]:.4f}, 損失={nonconvex_loss(z_grid[global_i]):.4f}")
plt.plot(z_grid, nonconvex_loss(z_grid), label="非凸な工程損失")
plt.scatter(nc_df["到達点"], nc_df["局所目的値"], color="tab:red", zorder=3, label="各初期値からの到達点")
plt.title("非凸問題では初期値により到達解が変わる")
plt.xlabel("標準化した工程条件 z"); plt.ylabel("損失指数")
plt.grid(True, alpha=0.3); plt.legend(); plt.tight_layout(); plt.show()

	初期値	到達点	局所目的値	成功
0	-1.6	-1.02178	-0.18198	True
1	-0.5	-1.02178	-0.18198	True
2	0.2	0.97669	0.17793	True
3	0.8	0.97669	0.17793	True
4	1.6	0.97669	0.17793	True

グリッドで確認した大域最良候補: z=-1.0224, 損失=-0.1820

結果の読み取り

ソルバーは各開始点で「成功」しても、左右の異なる谷へ到達し、目的値も異なります。非凸性が疑われる場合は、多点初期化、全域探索、混合整数最適化、シミュレーション、専門家による候補領域分割などを検討します。凸近似を使う場合は、元の制約に対する実行可能性と最適性ギャップを別に評価します。

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

1. まず凸性を確認する：実行可能領域と目的関数の構造を確認すると、解法だけでなく説明可能性が変わります。
2. 最適値より検算可能性を重視する：勾配、制約余裕、双対ギャップ、異なる初期値からの結果を残します。
3. 双対変数を投資判断へつなぐ：活性制約の緩和価値を円換算し、能力増強の候補を比較します。
4. 正則化を運用コストとして設計する：設定変更、センサー保守、モデル複雑度を目的へ明示的に組み込みます。
5. 凸でない現実を隠さない：段取り、整数判断、複数運転帯、不連続な品質判定があれば別のモデルや解法が必要です。

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

• 目的関数の合意：不良、再加工、電力、停止時間、納期遅延を共通単位へ換算し、部門間で承認する
• 制約の台帳化：設備上限、品質規格、安全余裕、変更可能幅について、根拠・責任者・更新日を管理する
• データ品質の確認：校正履歴、欠測、時刻同期、ロット差、条件変更の操作範囲を監査する
• オフライン検証：過去データでの再現、感度分析、外挿検知、異常入力時のフェイルセーフを試験する
• 現場での確認試験：推奨条件を段階的に試し、制御ばらつきを含めて規格余裕を評価する
• 運用監視：目的値、制約余裕、入力分布、予測誤差、再最適化理由を記録し、承認付きで反映する

まとめ

凸最適化は、単に速く解ける手法ではありません。安全な条件領域、損失の形、最適性条件、制約の価値を一貫して説明できる意思決定の枠組みです。本記事では、凸集合・凸関数から強双対性、制約付き最適化、Lasso、非凸問題の注意点までを、成形条件設計の流れで確認しました。実務では、モデルの凸性と現実の妥当性を分けて検証し、最適解を安全に実装できる条件へ翻訳することが重要です。

法人向けのご相談

数理工房では、製造条件最適化、生産計画、設備能力評価、品質予測と最適化の統合について、課題整理からPoC、本番運用設計まで支援しています。「目的関数をどう定義すればよいか」「制約や現場ルールをどうモデル化すればよいか」という段階からご相談いただけます。

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