"""Key metrics: torque, efficiency, backlash, torsional stiffness. Run: python key_metrics.py (reuses loaders/detectors from analysis.py) Writes km0_summary.png ... km4_backlash.png and key_metrics.txt next to this file. Angles are at the OUTPUT in degrees: motor travel (output-equivalent) = input turns * 360 / 30 = input turns * 12 output motion = MT6701 turns * 360 / 12 = MT6701 turns * 30 gearbox twist = motor travel - output motion """ import os import numpy as np import pandas as pd import matplotlib.pyplot as plt from matplotlib.lines import Line2D import analysis as A from analysis import (HERE, OUT, DPI, ARM_M, RATIO, SURF, INK, INK2, GRID, TILE, C_CH1, C_CH2, C_IN, C_SKIP, C_ALT) C_HEALTHY, C_DAMAGED = C_CH1, C_ALT def stops(k, d, c, s): """Per held contact: loading + unloading branches in (twist deg, torque N*m). Holds that contain a tooth skip are excluded (they are not elastic).""" t, inp, mt, T = d.time_s.values, d.in_turns.values, d.mt_turns.values, d.T_nm.values out = [] for _, r in c.iterrows(): if r.cmd == 0 or (k == "run1" and r.t1 > 41.3): # run1: only before skipping continue if ((s.t >= r.t0 - 0.1) & (s.t <= r.t1 + 0.3)).any(): continue a, b = np.searchsorted(t, r.t0), np.searchsorted(t, r.t1) sg = 1 if r.ch == "ch1" else -1 j = b # reversal = last stationary sample before the force drops while j > 0 and abs(inp[j] - inp[j - 1]) > 1e-3: j -= 1 ref = j seg = slice(max(a - 8, 0), b + 25) twist = sg * ((inp[seg] - inp[ref]) * 12 - (mt[seg] - mt[ref]) * 30) motor = sg * (inp[seg] - inp[ref]) * 12 outp = sg * (mt[seg] - mt[ref]) * 30 tq = sg * T[seg] idx = np.arange(seg.start, seg.stop) load = (idx >= a - 8) & (idx < a + np.argmax(sg * T[a:a + 30])) unload = idx >= j hold = r.F_hold * ARM_M fits = {} for name, m in (("load", load), ("unload", unload)): mm = m & (tq > 0.1 * hold) & (tq < 0.9 * hold) fits[name] = np.polyfit(twist[mm], tq[mm], 1)[0] if mm.sum() >= 4 else np.nan # output-motion dead band after the cell unloads iu = j while iu < len(T) - 1 and abs(T[iu]) > 0.015: iu += 1 # clean unload curve: twist measured back from the zero-torque point u = np.arange(j, iu + 1) u_tw = sg * ((inp[u] - inp[iu]) * 12 - (mt[u] - mt[iu]) * 30) u_tq = sg * T[u] mm = (u_tq > 0.1 * hold) & (u_tq < 0.9 * hold) k_u = np.polyfit(u_tw[mm], u_tq[mm], 1)[0] if mm.sum() >= 3 else np.nan out.append(dict(run=k, ch=r.ch, t0=r.t0, hold=hold, iq=r.iq, cmd=r.cmd, u_tw=u_tw, u_tq=u_tq, k_u=k_u, twist=twist, tq=tq, motor=motor, outp=outp, load=load, unload=unload, k_load=fits["load"], k_unload=fits["unload"], dep_motor=sg * (inp[j:j + 40] - inp[iu]) * 12 * -1, dep_out=sg * (mt[j:j + 40] - mt[iu]) * 30 * -1, dep_tq=np.abs(T[j:j + 40]), sample_deg=abs(inp[iu + 3] - inp[iu + 2]) * 12)) return out def main(): D = {k: A.load(k) for k in A.RUNS} C = {k: A.contacts(D[k]) for k in A.RUNS} S = {k: A.skips(D[k]) for k in ("run1", "run2")} cb = C["baseline"] kt = np.sum(cb.F_hold * ARM_M * cb.iq) / np.sum(cb.iq ** 2) ideal = kt * RATIO ST = stops("run1", D["run1"], C["run1"], S["run1"]) + stops("run2", D["run2"], C["run2"], S["run2"]) lines = [f"motor Kt {kt:.3f} N*m/A, ideal output {ideal:.2f} N*m/A"] # ---------------- KM1: torque timeline ------------------------------------- fig, axs = plt.subplots(1, 2, figsize=(14, 6), sharey=True, gridspec_kw=dict(width_ratios=[1.3, 1])) for ax, k, title in ((axs[0], "run1", "Run 1 (18:08), healthy at start"), (axs[1], "run2", "Run 2 (20:58), same gearbox, damaged")): c, s = C[k], S[k] for _, r in c[c.cmd > 0].iterrows(): sk = s[(s.t >= r.t0) & (s.t <= r.t1)] if len(sk): continue ax.scatter(r.t0, r.F_peak * ARM_M, s=110, marker="o", color=C_CH1 if r.ch == "ch1" else C_CH2, zorder=4) for _, e in s.iterrows(): ax.scatter(e.t, e.T_peak, s=45, marker="x", linewidths=2, color=C_CH1 if e.dir == "ch1" else C_CH2, zorder=5) ax.set_title(title, fontsize=13) ax.set_xlabel("time [s]") axs[0].axhline(1.47, color=INK2, ls="--", lw=1) axs[0].text(20.5, 1.50, "held 1.47 N·m without skipping", color=INK2, fontsize=11) axs[0].set_ylabel("peak output torque [N·m]") axs[0].set_ylim(0, 1.7) axs[1].legend(handles=[ Line2D([], [], marker="o", ls="", color=INK2, ms=10, label="held (no skip)"), Line2D([], [], marker="x", ls="", color=INK2, ms=9, mew=2, label="tooth skip"), Line2D([], [], marker="s", ls="", color=C_CH1, ms=10, label="CH1 direction"), Line2D([], [], marker="s", ls="", color=C_CH2, ms=10, label="CH2 direction")], loc="upper right") fig.suptitle("TORQUE: ratchets at ~1.5 N·m healthy, then ~1.1, then 0.3–0.7 N·m once damaged", x=0.01, ha="left", fontsize=15, fontweight="bold", color=INK) fig.tight_layout() fig.savefig(os.path.join(OUT, "km1_torque.png"), dpi=DPI) # ---------------- KM2: efficiency ------------------------------------------ fig, ax = plt.subplots(figsize=(11, 6.5)) rows = [] for k in ("run1", "run2"): c = C[k] c = c[(c.cmd > 0) & (c.iq > 0.05)] if k == "run1": c = c[c.t1 < 41.3] for _, r in c.iterrows(): eff = r.F_hold * ARM_M / (ideal * r.iq) * 100 rows.append(dict(run=k, ch=r.ch, Tq=r.F_hold * ARM_M, iq=r.iq, eff=eff)) E = pd.DataFrame(rows) style = {("run1", "ch1"): (C_CH1, "o"), ("run1", "ch2"): (C_CH2, "o"), ("run2", "ch1"): (C_CH1, "D"), ("run2", "ch2"): (C_CH2, "D")} for (k, ch), g in E.groupby(["run", "ch"]): col, mk = style[(k, ch)] face = col if k == "run1" else "none" ax.scatter(g.Tq, g.eff, s=130, marker=mk, facecolors=face, edgecolors=col, linewidths=2.2, zorder=4) for n, (_, r) in enumerate(g.iterrows()): dy = -16 if (k == "run1" and ch == "ch1" and r.Tq < 1.3) else -4 ax.annotate(f"{r.iq:.2f} A", (r.Tq, r.eff), textcoords="offset points", xytext=(9, dy), color=INK2, fontsize=10) ax.set_xlabel("held output torque [N·m]") ax.set_ylabel("static efficiency [%]\n= held torque / (30 × motor torque)") ax.set_xlim(0, 1.5) ax.set_ylim(0, 75) ax.legend(handles=[ Line2D([], [], marker="o", ls="", color=C_CH1, ms=10, label="run 1 healthy, CH1 dir"), Line2D([], [], marker="o", ls="", color=C_CH2, ms=10, label="run 1 healthy, CH2 dir"), Line2D([], [], marker="D", ls="", mfc="none", mec=C_CH1, mew=2, ms=9, label="run 2 damaged, CH1 dir"), Line2D([], [], marker="D", ls="", mfc="none", mec=C_CH2, mew=2, ms=9, label="run 2 damaged, CH2 dir")], loc="upper right") ax.set_title("EFFICIENCY: 43–55 % healthy (stalled), the worn direction drops to 8–20 %") fig.tight_layout() fig.savefig(os.path.join(OUT, "km2_efficiency.png"), dpi=DPI) lines.append("\nEfficiency per hold:\n" + E.round(3).to_string()) # ---------------- KM3: stiffness ------------------------------------------- # Unloading branch only: starts from rest (0.1 s speed slew), so it is quasi-static; # the loading branch carries the impact of hitting the stop at speed. fig, ax = plt.subplots(figsize=(11, 7)) ks = {"run1": [], "run2": []} for st in ST: col = C_HEALTHY if st["run"] == "run1" else C_DAMAGED mk = "o" if st["ch"] == "ch1" else "s" ax.plot(st["u_tw"], st["u_tq"], color=col, lw=1.8, marker=mk, ms=4.5, alpha=0.9) if np.isfinite(st["k_u"]) and st["k_u"] > 0: ks[st["run"]].append(dict(ch=st["ch"], hold=st["hold"], k=st["k_u"])) for k, col, xmax in (("run1", C_HEALTHY, 2.3), ("run2", C_DAMAGED, 2.4)): K = pd.DataFrame(ks[k]) kd, lo, hi = K.k.median(), K.k.min(), K.k.max() xx = np.linspace(0, xmax, 10) ax.plot(xx, kd * xx, color=col, ls="--", lw=1.3, alpha=0.8) lines.append(f"stiffness {k} (unload): median {kd:.3f} N*m/deg = {kd*57.3:.1f} N*m/rad, " f"range {lo*57.3:.1f}-{hi*57.3:.1f}, n={len(K)}") ks[k + "_med"] = kd k1, k2 = ks["run1_med"], ks["run2_med"] ax.text(2.75, 1.12, f"healthy: {k1:.2f} N·m/°\n= {k1 * 57.3:.0f} N·m/rad", color=C_HEALTHY, fontsize=13, fontweight="bold") ax.text(0.9, 0.06, f"damaged: {k2:.2f} N·m/° = {k2 * 57.3:.0f} N·m/rad", color=C_DAMAGED, fontsize=13, fontweight="bold") ax.annotate("damaged teeth peak at 0.5–0.7 N·m, then give way:\n" "the holds were parked part-way over the tooth tips\n(read right → left as the motor backs off)", xy=(2.65, 0.70), xytext=(2.9, 0.78), color=INK2, fontsize=11, arrowprops=dict(arrowstyle="->", color=INK2)) ax.set_xlim(0, 5) ax.set_ylim(0, 1.5) ax.set_xlabel("gearbox twist [deg at the output] (motor/30 minus measured arm motion)") ax.set_ylabel("output torque [N·m]") ax.legend(handles=[Line2D([], [], color=C_HEALTHY, lw=3, label="run 1 healthy"), Line2D([], [], color=C_DAMAGED, lw=3, label="run 2 damaged"), Line2D([], [], color=INK2, marker="o", ls="", label="CH1 direction"), Line2D([], [], color=INK2, marker="s", ls="", label="CH2 direction")], loc="upper left") ax.set_title(f"TORSIONAL STIFFNESS: halved after skipping ({k1 * 57.3:.0f} → {k2 * 57.3:.0f} N·m/rad)") fig.tight_layout() fig.savefig(os.path.join(OUT, "km3_stiffness.png"), dpi=DPI) # ---------------- KM4: backlash -------------------------------------------- fig, axs = plt.subplots(2, 1, figsize=(11, 8), sharex=True, gridspec_kw=dict(height_ratios=[1, 1.3])) res = [] for st in ST: col = C_HEALTHY if st["run"] == "run1" else C_DAMAGED x = st["dep_motor"] axs[0].plot(x, st["dep_tq"], color=col, lw=1.6, alpha=0.85, marker="o", ms=3) axs[1].plot(x, st["dep_out"], color=col, lw=1.6, alpha=0.85, marker="o", ms=3) res.append(st["sample_deg"]) xx = np.linspace(0, 3, 10) axs[1].plot(xx, xx, color=INK2, ls="--", lw=1.2) axs[1].text(1.6, 1.15, "arm follows motor 1:1\n(zero backlash)", color=INK2, fontsize=11) axs[1].plot([0, 1.0, 3], [0, 0, 2.0], color=C_SKIP, ls=":", lw=1.6) axs[1].text(1.05, -0.35, "what 1° of backlash would look like", color=C_SKIP, fontsize=11) for ax in axs: ax.axvline(0, color=INK, lw=1) ax.axvspan(-5, 0, color=GRID, alpha=0.35, lw=0) axs[0].text(-4.8, axs[0].get_ylim()[1] * 0.85, "motor reversing,\ngearbox unwinding (elastic)", color=INK2) axs[0].text(0.08, axs[0].get_ylim()[1] * 0.85, "load cell reads 0", color=INK) axs[0].set_ylabel("torque on stop [N·m]") axs[1].set_ylabel("arm motion off the stop [deg]") axs[1].set_xlabel("motor travel since the cell unloaded [deg, output-equivalent]") axs[1].set_xlim(-5, 3) axs[1].set_ylim(-0.6, 3) axs[1].legend(handles=[Line2D([], [], color=C_HEALTHY, lw=3, label="run 1 healthy (300 RPM)"), Line2D([], [], color=C_DAMAGED, lw=3, label="run 2 damaged (150 RPM)")], loc="upper left") axs[0].set_title(f"BACKLASH: arm moves the instant the load is gone; dead band < {max(res):.2f}° (1 sample)") fig.tight_layout() fig.savefig(os.path.join(OUT, "km4_backlash.png"), dpi=DPI) # ---------------- KM5: hysteresis ------------------------------------------ # Per held stop: load (twist up) -> hold (motor stalled, torque relaxes) -> unload. # The arm sits on the same cell at both zero-torque ends, so the MT6701 scale # error cancels and the twist offset there is the residual set. def mono_interp(lvl, tq, tw): tq_m = np.maximum.accumulate(tq) return np.interp(lvl, tq_m, tw) HY = [] for k in ("run1", "run2"): d, c, s = D[k], C[k], S[k] t, inp, mt, T = d.time_s.values, d.in_turns.values, d.mt_turns.values, d.T_nm.values for _, r in c.iterrows(): if r.cmd == 0 or (k == "run1" and r.t1 > 41.3): continue if ((s.t >= r.t0 - 0.1) & (s.t <= r.t1 + 0.3)).any(): continue sg = 1 if r.ch == "ch1" else -1 a, b = np.searchsorted(t, r.t0), np.searchsorted(t, r.t1) tq, tw = sg * T, sg * (inp * 12 - mt * 30) i0 = a while i0 > 0 and tq[i0] > 0.03: i0 -= 1 ipk = a + np.argmax(tq[a:a + 40]) j = b while j > 0 and abs(inp[j] - inp[j - 1]) > 1e-3: j -= 1 i1 = b while tq[i1] > 0.03: i1 += 1 hold = r.F_hold * ARM_M lo_lvl = 0.06 ref = mono_interp(lo_lvl, tq[i0:ipk + 1], tw[i0:ipk + 1]) L = lambda lvl: mono_interp(lvl, tq[i0:ipk + 1], tw[i0:ipk + 1]) - ref U = lambda lvl: mono_interp(lvl, tq[j:i1 + 1][::-1], tw[j:i1 + 1][::-1]) - ref path = np.arange(i0, i1 + 1) x, y = tw[path] - ref, tq[path] loop = np.trapezoid(y, np.radians(x)) # dissipated per cycle [J] w_in = np.trapezoid(tq[i0:j + 1], np.radians(tw[i0:j + 1])) # work put in HY.append(dict(run=k, ch=r.ch, hold=hold, x=x, y=y, n_load=ipk - i0, n_hold=j - i0, set=U(lo_lvl) - L(lo_lvl), w50=U(0.5 * hold) - L(0.5 * hold), loss_mJ=loop * 1e3, loss_pct=100 * loop / w_in, relax=tq[ipk] - tq[j], relax_pct=100 * (tq[ipk] - tq[j]) / tq[ipk], hold_s=t[j] - t[ipk])) H = pd.DataFrame([{kk: v for kk, v in h.items() if kk not in ("x", "y")} for h in HY]) lines.append("\nHysteresis per stop (deg at output; w50 = load-vs-unload twist gap at half torque, " "set = residual twist at ~0 torque):\n" + H.round(3).to_string()) hh = H[H.run == "run1"] hd1 = H[(H.run == "run2") & (H.ch == "ch1")] hd2 = H[(H.run == "run2") & (H.ch == "ch2")] lines.append(f"healthy: w50 median {hh.w50.median():.3f} deg, set median {hh.set.median():.3f} deg, " f"loss {hh.loss_pct.median():.0f}% ({hh.loss_mJ.median():.1f} mJ), relax {hh.relax_pct.median():.0f}% in ~2 s") fig, axs = plt.subplots(1, 2, figsize=(15, 6.8), sharey=True) for h in HY: ax = axs[0] if h["run"] == "run1" else axs[1] col = C_CH1 if h["ch"] == "ch1" else C_CH2 ax.plot(h["x"], h["y"], color=col, lw=1.8, marker="o", ms=3, alpha=0.9) # annotate one healthy loop's phases h0 = max([h for h in HY if h["run"] == "run1"], key=lambda h: h["hold"]) xl, yl = h0["x"], h0["y"] axs[0].annotate("load ↗", xy=(xl[h0["n_load"] // 2], yl[h0["n_load"] // 2]), xytext=(0.15, 1.05), color=INK, fontsize=12, arrowprops=dict(arrowstyle="->", color=INK2)) axs[0].annotate("hold: motor stalled,\ntorque relaxes ↓ (creep)", xy=(xl[h0["n_hold"]], yl[h0["n_hold"]]), xytext=(2.45, 0.95), color=INK, fontsize=12, arrowprops=dict(arrowstyle="->", color=INK2)) iu = h0["n_hold"] + int(np.argmax(yl[h0["n_hold"]:] < 0.5 * h0["hold"])) # half-torque on unload axs[0].annotate("unload ↙", xy=(xl[iu], yl[iu]), xytext=(2.2, 0.4), color=INK, fontsize=12, arrowprops=dict(arrowstyle="->", color=INK2)) axs[0].set_title(f"Run 1, healthy: narrow loop that closes\n" f"width {hh.w50.median():.2f}° at half torque · set {hh.set.median():+.2f}° · " f"{hh.loss_pct.median():.0f}% work lost", fontsize=12.5) axs[1].set_title(f"Run 2, damaged: fat CH2 loop (teeth over the tips)\n" f"work lost: CH1 {hd1.loss_pct.median():.0f}% · CH2 {hd2.loss_pct.median():.0f}% · " f"CH2 relaxed {hd2.relax_pct.median():.0f}% while held", fontsize=12.5) axs[0].set_xlim(-0.2, 3.6) axs[1].set_xlim(-0.3, 5.0) for ax in axs: ax.set_xlabel("gearbox twist from start of loading [deg at the output]") axs[0].set_ylabel("output torque on stop [N·m]") axs[0].set_ylim(-0.03, 1.55) axs[1].legend(handles=[Line2D([], [], color=C_CH1, lw=3, label="CH1 direction"), Line2D([], [], color=C_CH2, lw=3, label="CH2 direction")], loc="upper right") fig.suptitle("HYSTERESIS: torque–twist loops at each stop (arm/cell motion subtracted)", x=0.01, ha="left", fontsize=15, fontweight="bold", color=INK) fig.tight_layout() fig.savefig(os.path.join(OUT, "km5_hysteresis.png"), dpi=DPI) # ---------------- KM0: summary card ---------------------------------------- k1, k2 = ks["run1_med"] * 57.3, ks["run2_med"] * 57.3 e1 = E[E.run == "run1"].eff tiles = [ ("TORQUE", "1.47 N·m", "held before first skip\n(ratchets at 1.48)", "damaged: 0.3–0.7 N·m", C_CH2), ("EFFICIENCY", f"{e1.min():.0f}–{e1.max():.0f} %", "static, stalled\n(falls as load rises)", "worn side: 8–20 %", C_IN), ("BACKLASH", f"< {max(res):.1f}°", "below 1-sample resolution\n(≈ 0 at the gears)", "slow test needed for exact value", C_ALT), ("STIFFNESS", f"{k1:.0f} N·m/rad", f"= {k1 / 57.3:.2f} N·m per degree\n(stiffer at higher load)", f"damaged: {k2:.0f} N·m/rad", C_CH1), ("HYSTERESIS", f"{hh.loss_pct.median():.0f} % loss", f"loop width {hh.w50.median():.2f}° at half torque\nset {hh.set.median():.1f}° · relaxes {hh.relax_pct.median():.0f}% in 2 s", f"damaged CH2: {hd2.loss_pct.median():.0f}% loss", C_CH2), ] fig, axs = plt.subplots(1, 5, figsize=(20, 4.6)) for ax, (name, big, sub, note, col) in zip(axs, tiles): ax.axis("off") ax.add_patch(plt.Rectangle((0, 0), 1, 1, transform=ax.transAxes, facecolor=TILE, edgecolor=GRID, lw=1)) ax.add_patch(plt.Rectangle((0, 0.96), 1, 0.04, transform=ax.transAxes, facecolor=col, lw=0)) ax.text(0.07, 0.82, name, transform=ax.transAxes, color=INK2, fontsize=13, fontweight="bold") ax.text(0.07, 0.55, big, transform=ax.transAxes, color=INK, fontsize=30, fontweight="bold") ax.text(0.07, 0.30, sub, transform=ax.transAxes, color=INK2, fontsize=11.5, va="center") ax.text(0.07, 0.08, note, transform=ax.transAxes, color=C_SKIP, fontsize=11.5) fig.suptitle("Harmonic Pancake RevB 30:1 (PETG flex / PLA): healthy numbers, run 1", x=0.01, ha="left", fontsize=15, fontweight="bold", color=INK) fig.tight_layout() fig.savefig(os.path.join(OUT, "km0_summary.png"), dpi=DPI) txt = "\n".join(lines) open(os.path.join(HERE, "key_metrics.txt"), "w").write(txt) print(txt) if __name__ == "__main__": main()