はじめに

2016年、エアバスのエンジニアは、AIが提示した「骨のような有機構造」の仕切り板に絶句しました。従来比 45%軽く、強度は同等以上。これが ジェネレーティブ・デザイン（Generative Design） が世界に投げかけた最初の衝撃です。設計者は「描く人」から「AIが提案した数千案から選ぶキュレーター」になりつつあります。

「でもAutodesk Fusion / nTopologyみたいなツールがないと無理でしょ？」――いえ、Python + numpy/scipy + pymoo + trimesh などのOSSを組み合わせれば、ジェネレーティブデザインの コアロジック（トポロジー最適化＋多目的最適化＋メッシュ出力） を手元のPCで一通り体験できます。

本記事では、note記事のテーマを開発者目線で落とし込み、Pythonでジェネレーティブデザイン最小ひな形を構築する具体的な手順 を解説します。

参考にした解説記事：「人間には絶対に思いつかない形」——ジェネレーティブ・デザインが、製造業の"設計"を再発明する | AI Robotics Quantum Lab

ゴールは以下のとおりです。

2D片持ち梁（MBB beam）に対し SIMP法でトポロジー最適化 を解く
制約条件（体積率・荷重位置）を変えて 複数案を一気に生成 する
pymooのNSGA-II で「軽さ ↔ 剛性 ↔ 製造容易性」を多目的最適化する
結果を等高線抽出して STL／PNG で書き出す（3Dプリント連携）
機械学習サロゲート、3D拡張、CAD/3Dプリンタ連携への道筋を整理する

対象読者は、Pythonと数値計算の基本に触れたことがある、機械設計・CAE・3Dプリンタに興味のあるエンジニアです。

1. 全体像 ― 何を作るのか

[設計領域 + 荷重・拘束条件]
       │
       ▼  SIMP法トポロジー最適化（FEM + 密度設計変数）
[最適密度マップ ρ(x,y)]
       │
       ├──► 単目的: 体積制約下のコンプライアンス最小化
       │
       └──► 多目的: 重量↓ + 剛性↑ + 製造性(オーバーハング小) ↓
              │
              ▼  pymoo NSGA-II
       [パレートフロント上の候補群]
              │
              ▼  Marching Squares → 等高線 → STL押し出し
       [製造可能な3Dモデル群]
              │
              ▼
       matplotlib / PyVista で可視化

「問いを与え、AIが何千案も出す。人間は選ぶ」を、コードでそのまま再現します。

2. 環境構築

2.1 Python仮想環境

mkdir generative-design && cd generative-design
python3 -m venv .venv
source .venv/bin/activate
pip install --upgrade pip

2.2 必要ライブラリ

# 数値計算
pip install numpy scipy matplotlib

# 最適化
pip install pymoo                  # 多目的（NSGA-II）

# 形状処理・3D
pip install scikit-image           # 等高線抽出
pip install trimesh                # STL生成
pip install pyvista                # 3D可視化

# UI
pip install streamlit

2.3 プロジェクト構成

generative-design/
├── src/
│   ├── topopt.py          # SIMP法本体
│   ├── multi_obj.py       # NSGA-II ラッパー
│   ├── mesh_export.py     # STL書き出し
│   └── viz.py
├── outputs/
└── app.py                 # Streamlit

3. SIMP法によるトポロジー最適化（最小実装）

SIMP（Solid Isotropic Material with Penalization）は、各セルに密度 ρ∈[0,1] を割り当て、ρ^p で剛性を罰則化することで「材料があるか／ないか」をなめらかに最適化する手法です。

# src/topopt.py
import numpy as np
from scipy.sparse import coo_matrix
from scipy.sparse.linalg import spsolve


def topopt_2d(nelx: int = 60, nely: int = 30,
              volfrac: float = 0.4, penal: float = 3.0,
              rmin: float = 1.5, max_iter: int = 50,
              load_pos: tuple[int, int] | None = None):
    """2D片持ち梁のSIMPトポロジー最適化。
    Returns: 密度マップ ρ shape=(nely, nelx)
    """
    # 要素剛性マトリクス（4節点平面応力、E=1, nu=0.3）
    KE = _element_stiffness()
    edofMat, iK, jK = _build_dof(nelx, nely)

    x = np.full((nely, nelx), volfrac, dtype=float)

    for it in range(max_iter):
        # 全体剛性マトリクスを組み立て
        sK = ((KE.flatten()[np.newaxis]).T * (x.T.flatten() ** penal)).flatten(order="F")
        K = coo_matrix((sK, (iK, jK)), shape=((nelx + 1) * (nely + 1) * 2,) * 2).tocsr()

        # 境界条件
        ndof = (nelx + 1) * (nely + 1) * 2
        fixed = np.union1d(np.arange(0, 2 * (nely + 1), 2),     # 左端x固定
                            np.arange(1, 2 * (nely + 1), 2))   # 左端y固定
        free = np.setdiff1d(np.arange(ndof), fixed)

        # 荷重
        F = np.zeros(ndof)
        if load_pos is None:
            load_pos = (nelx, nely // 2)   # 右下/右中央に下向き荷重
        node = load_pos[0] * (nely + 1) + load_pos[1]
        F[2 * node + 1] = -1.0

        U = np.zeros(ndof)
        U[free] = spsolve(K[free, :][:, free], F[free])

        # 要素ごとの感度
        ce = np.einsum("ij,jk,ki->i", U[edofMat], KE, U[edofMat].T)
        ce = ce.reshape((nelx, nely)).T

        compliance = float((x ** penal * ce).sum())
        dc = -penal * (x ** (penal - 1)) * ce

        # 感度フィルタ（製造性のための平滑化）
        dc = _sensitivity_filter(x, dc, rmin)

        # OC法で更新
        x = _oc_update(x, dc, volfrac)

        if it % 10 == 0:
            print(f"iter={it:3d}  compliance={compliance:.3f}  vol={x.mean():.3f}")

    return x


# --- 下請け関数たち ---
def _element_stiffness():
    E, nu = 1.0, 0.3
    k = np.array([0.5 - nu/6, 0.125 + nu/8, -0.25 - nu/12, -0.125 + 3*nu/8,
                  -0.25 + nu/12, -0.125 - nu/8, nu/6, 0.125 - 3*nu/8])
    KE = (E / (1 - nu**2)) * np.array([
        [k[0], k[1], k[2], k[3], k[4], k[5], k[6], k[7]],
        [k[1], k[0], k[7], k[6], k[5], k[4], k[3], k[2]],
        [k[2], k[7], k[0], k[5], k[6], k[3], k[4], k[1]],
        [k[3], k[6], k[5], k[0], k[7], k[2], k[1], k[4]],
        [k[4], k[5], k[6], k[7], k[0], k[1], k[2], k[3]],
        [k[5], k[4], k[3], k[2], k[1], k[0], k[7], k[6]],
        [k[6], k[3], k[4], k[1], k[2], k[7], k[0], k[5]],
        [k[7], k[2], k[1], k[4], k[3], k[6], k[5], k[0]],
    ])
    return KE


def _build_dof(nelx, nely):
    nodenrs = np.arange((nelx + 1) * (nely + 1)).reshape((nelx + 1, nely + 1)).T
    edofVec = (2 * nodenrs[:-1, :-1] + 2).flatten(order="F")
    edofMat = np.tile(edofVec, (8, 1)).T + np.array(
        [0, 1, 2 * nely + 2, 2 * nely + 3, 2 * nely, 2 * nely + 1, -2, -1]
    )
    iK = np.kron(edofMat, np.ones((8, 1))).flatten()
    jK = np.kron(edofMat, np.ones((1, 8))).flatten()
    return edofMat, iK, jK


def _sensitivity_filter(x, dc, rmin):
    nely, nelx = x.shape
    H = np.zeros_like(x)
    for i in range(nelx):
        for j in range(nely):
            i1, i2 = max(0, i - int(np.ceil(rmin))), min(nelx, i + int(np.ceil(rmin)) + 1)
            j1, j2 = max(0, j - int(np.ceil(rmin))), min(nely, j + int(np.ceil(rmin)) + 1)
            ii, jj = np.meshgrid(np.arange(i1, i2), np.arange(j1, j2), indexing="ij")
            w = np.maximum(0, rmin - np.sqrt((ii - i) ** 2 + (jj - j) ** 2))
            H[j, i] = (w * x[j1:j2, i1:i2].T * dc[j1:j2, i1:i2].T).sum() / (x[j, i] * w.sum() + 1e-9)
    return H


def _oc_update(x, dc, volfrac, move: float = 0.2):
    l1, l2 = 0.0, 1e9
    while l2 - l1 > 1e-4:
        lmid = 0.5 * (l1 + l2)
        x_new = np.clip(
            np.maximum(1e-3, np.maximum(x - move, np.minimum(1.0, np.minimum(x + move,
                       x * np.sqrt(-dc / lmid + 1e-9))))),
            0.001, 1.0,
        )
        if x_new.mean() > volfrac:
            l1 = lmid
        else:
            l2 = lmid
    return x_new

💡 この実装は 教育用に簡略化 した形です。本格的にはOlivier Sigmund の "Top88" 系コードや、topopt/SciPy + cvxopt ベースのライブラリを使うと素直です。

実行：

import matplotlib.pyplot as plt
from src.topopt import topopt_2d

rho = topopt_2d(nelx=80, nely=40, volfrac=0.4, max_iter=50)
plt.imshow(1 - rho, cmap="gray")
plt.title("Generative shape (SIMP)")
plt.colorbar()
plt.savefig("outputs/simp_result.png", dpi=150)

「片持ち梁の有機トラス」が黒い線として現れれば成功です。「人間が描かない、骨のような形」 が、コード50行ちょっとから出てきます。

4. 「数千案を生み出す」ジェネレーティブ・スイープ

ジェネレーティブ・デザインの本質は 「条件を変えた最適化を大量に回す」 こと。volfrac・荷重位置・フィルタ半径を変えて一気に走らせます。

import itertools, os
from src.topopt import topopt_2d

os.makedirs("outputs/sweep", exist_ok=True)

cases = list(itertools.product(
    [0.25, 0.35, 0.45, 0.55],         # 体積率
    [(80, 20), (80, 30), (80, 40)],   # 荷重位置
    [1.2, 1.5, 2.0],                  # フィルタ半径
))

for i, (vf, lp, rmin) in enumerate(cases):
    rho = topopt_2d(nelx=80, nely=40, volfrac=vf,
                    load_pos=lp, rmin=rmin, max_iter=40)
    plt.imsave(f"outputs/sweep/{i:03d}_vf{vf}_lp{lp}_r{rmin}.png",
               1 - rho, cmap="gray")

得られた数十枚を眺めれば、「条件を変えるだけで全く違う構造が立ち上がる」 という、ジェネレーティブの醍醐味が見えます。

5. 多目的最適化（pymoo NSGA-II）

「軽さ ↔ 剛性 ↔ 製造性」のように、トレードオフの異なる目的を同時に解く のがジェネレーティブ・デザインの真骨頂です。pymoo の NSGA-II を使えば、パレートフロント上の候補群 が得られます。

# src/multi_obj.py
import numpy as np
from pymoo.core.problem import ElementwiseProblem
from pymoo.algorithms.moo.nsga2 import NSGA2
from pymoo.operators.crossover.sbx import SBX
from pymoo.operators.mutation.pm import PM
from pymoo.operators.sampling.rnd import FloatRandomSampling
from pymoo.optimize import minimize

from src.topopt import topopt_2d


class BracketDesign(ElementwiseProblem):
    """設計変数: [volfrac, load_y(0~1), rmin]
       目的1: 体積使用量（小さいほど良）
       目的2: コンプライアンス（小さいほど剛性高）
       目的3: 製造性ペナルティ（中央集中度合いを近似）
    """

    def __init__(self):
        super().__init__(n_var=3, n_obj=3,
                         xl=np.array([0.20, 0.10, 1.0]),
                         xu=np.array([0.60, 0.90, 2.5]))

    def _evaluate(self, x, out, *args, **kwargs):
        volfrac, ly, rmin = x
        nelx, nely = 60, 30
        load_y = int(ly * nely)
        rho = topopt_2d(nelx=nelx, nely=nely, volfrac=volfrac,
                        rmin=rmin, max_iter=20,
                        load_pos=(nelx, load_y))
        vol_used = float(rho.mean())
        # ピクセル単位の応力評価は省略し、密度の濃淡で疑似コンプライアンス
        compliance = float(((1 - rho) ** 2).sum())
        manufacturability = float(((rho > 0.5).astype(float).std()))
        out["F"] = [vol_used, compliance, manufacturability]


def run_nsga(pop: int = 12, n_gen: int = 5):
    problem = BracketDesign()
    algo = NSGA2(
        pop_size=pop,
        sampling=FloatRandomSampling(),
        crossover=SBX(prob=0.9, eta=15),
        mutation=PM(eta=20),
        eliminate_duplicates=True,
    )
    res = minimize(problem, algo, ("n_gen", n_gen),
                   seed=42, verbose=True)
    return res

res = run_nsga(pop=12, n_gen=5)
print("Pareto front:")
print(res.F)
print("Design variables:")
print(res.X)

⚠️ NSGA-II の評価1回ごとに SIMP を解くので、本格化する際は max_iter を絞る or サロゲートモデルを挟む ことが必須です（次節）。

6. 機械学習サロゲートモデル

NSGA-IIで毎回FEMを解くのは重いです。事前に少量のSIMP結果を集め、ニューラルネットで「条件→目的値」を学習 すれば、評価が数千倍速くなります。

# scikit-learnでの最小サンプル
import numpy as np
from sklearn.ensemble import GradientBoostingRegressor

# X: (volfrac, ly, rmin)  Y: (vol_used, compliance, manuf)
X = np.random.rand(200, 3)      # 実際はSIMPを200回回して集める
Y = np.random.rand(200, 3)

models = [GradientBoostingRegressor().fit(X, Y[:, j]) for j in range(3)]

def surrogate(x):
    return np.array([m.predict(x.reshape(1, -1))[0] for m in models])

NSGA-II の _evaluate を surrogate(x) に置き換えると、1案あたり数ms で評価できます。「ML surrogate + Active Learning でSIMPを呼ぶ回数を最小化」が現代のジェネレーティブ・デザインの常套パターンです。

7. 結果を3Dプリンタ向けにSTL出力

2D密度マップ → 等高線抽出 → 押し出してSTL化、の流れです。

# src/mesh_export.py
import numpy as np
import trimesh
from skimage import measure


def density_to_stl(rho: np.ndarray, thickness_mm: float = 5.0,
                   cell_mm: float = 1.0, out_path: str = "outputs/part.stl"):
    binary = (rho > 0.5).astype(np.uint8)
    contours = measure.find_contours(binary, 0.5)

    meshes = []
    for c in contours:
        pts2d = np.column_stack([c[:, 1] * cell_mm, c[:, 0] * cell_mm])
        if len(pts2d) < 3:
            continue
        # 単純押し出し（凸性を仮定）
        poly = trimesh.path.polygons.paths_to_polygons([pts2d])
        if not poly:
            continue
        mesh = trimesh.creation.extrude_polygon(poly[0], height=thickness_mm)
        meshes.append(mesh)

    combined = trimesh.util.concatenate(meshes)
    combined.export(out_path)
    print(f"saved: {out_path}")

from src.mesh_export import density_to_stl
density_to_stl(rho, thickness_mm=5.0, cell_mm=2.0)

outputs/part.stl を Cura / PrusaSlicer に投げれば、そのまま3Dプリンターで実物化 できます。家にFDMプリンタがあれば、AIが設計した「骨のような梁」を手のひらで握れる体験はかなり強烈です。

8. Streamlitで触れるダッシュボード

# app.py
import streamlit as st
import numpy as np
import matplotlib.pyplot as plt
from src.topopt import topopt_2d


st.set_page_config(layout="wide")
st.title("Generative Design – SIMP Playground")

c1, c2 = st.columns(2)
with c1:
    volfrac = st.slider("Volume fraction", 0.1, 0.7, 0.4, step=0.05)
    rmin = st.slider("Filter radius", 1.0, 3.0, 1.5, step=0.1)
    load_y_ratio = st.slider("Load position (y)", 0.0, 1.0, 0.5, step=0.05)
    max_iter = st.slider("Iterations", 10, 80, 40, step=5)

if st.button("Generate"):
    nelx, nely = 80, 40
    rho = topopt_2d(nelx, nely, volfrac=volfrac, rmin=rmin,
                    max_iter=max_iter,
                    load_pos=(nelx, int(load_y_ratio * nely)))
    fig, ax = plt.subplots()
    ax.imshow(1 - rho, cmap="gray")
    ax.set_axis_off()
    c2.pyplot(fig)
    c2.metric("Volume used", f"{rho.mean():.3f}")

streamlit run app.py

スライダーを動かすたびに「全く違う有機構造」が立ち上がる、「設計はパラメータの問いかけだ」 という感覚をその場で体験できるツールになります。

9. 3Dへの拡張・実機CAD連携

3Dトポロジー最適化：voxel格子＋3D-FEM。SfePy, FEniCS, MFEM などで本格的に。GPU化なら JAX-FEM も選択肢。
応力以外の物理：熱、流体、振動を加える「マルチフィジクス・ジェネレーティブ」。
CAD連携：Autodesk Fusion API、SolidWorks API、FreeCAD（Pythonで操作可）。
3Dプリンタ最適化：自重支持の限界（オーバーハング角）を制約に組み込んだ「自立設計」最適化。
マテリアルズ・インフォマティクス：素材物性（ヤング率・降伏応力）まで設計変数化し、構造と材料を 同時に 最適化する。

10. 製造業導入の運用上の注意

出力された有機形状は 金型成形が難しい ことが多いです。射出成形なら抜き勾配・パーティングライン制約を加え、3Dプリント前提なら印刷方位を制約に。
解析と現実のギャップ：FEMの簡略化（線形・等方）が破綻するケース（座屈・疲労・接合部）を必ず実機試験で詰める。
規制対応：航空・医療・建築では、設計プロセスのトレーサビリティ が法定要件。最適化の入力・結果・人間の判断ログを残す仕組みを設計初期から組み込む。
知財：AIが生成した形状の著作権・特許は国・分野で扱いが分かれます。設計データ管理を契約・社内規程で明確化。
人間の役割：最終判断は人間。AIの提案を盲信しないレビュー体制（DRレビュー、FMEA）を残す。

11. 拡張アイデア

拡散モデル × 構造設計：DDPMで「構造画像」を直接生成し、SIMP結果を教師に学習
強化学習：荷重・拘束をエージェントに与え、構造を逐次配置するポリシーを学習
微分可能シミュレーション：DiffTaichi / JAX-FEMでエンドツーエンド最適化
CAEクラウド連携：Ansys Cloud / Autodesk Forge と組み合わせた現業向けPoC
生体模倣ライブラリ：骨・葉脈・サンゴの自然構造パターンを「事前モデル」として組み込む

まとめ

ジェネレーティブ・デザインの核は 「トポロジー最適化 × 多目的最適化 × 大量試行」
本記事のひな形で、SIMP法 → 条件スイープ → NSGA-II多目的 → STL書き出し → Streamlit までを最短ルートで体験できる
ML surrogate + 3D拡張 + CAD連携で、実務に届くレベルまでステップアップ可能
設計者の本質は「線を引くこと」から「問いを立て、AIの提案から選ぶこと」に書き換わる。エンジニアはその新しい仕事を回す 道具立てとロジック を作れる立場にある

まずは片持ち梁のSIMPを走らせ、「人間が思いつかない骨のような形」がコードから立ち上がる瞬間を体感してみてください。製造業の"設計"を再発明する 流れに自分の手で乗っていく入口は、自分のノートPCの中にあります。