Port MATLAB geometry+process optimization to Python

docs/tutorial.md

@@ -131,7 +131,59 @@ See [`python/examples/`](../python/examples/) for a runnable script for each.
 
 ## Optimization
 
-Design optimization (resolution / noise PSD goal functions with parameter and nonlinear constraints) is implemented in the MATLAB version but has not yet been ported to Python. Track [GitHub Issues](https://github.com/MicrosystemsLab/PiezoD/issues) for progress, or use the MATLAB implementation when optimization is required.
+Design optimization is built around three pieces:
+
+- A **goal callable** that maps a cantilever to the scalar to minimize. PiezoD ships factories that match the MATLAB tutorial's units (pN, nm, pN/sqrt(Hz), etc.):
+
+```python
+from piezod import (
+    force_resolution_goal,
+    displacement_resolution_goal,
+    force_noise_density_goal,
+    surface_stress_resolution_goal,
+)
+```
+
+- Optional **parameter_constraints** that override the default state-variable bounds (e.g. clamp thickness, cap bias voltage). Keys are `min_<name>` / `max_<name>` for any state variable returned by `c.optimization_state_vars()`.
+
+- Optional **metric_constraints** -- inequality constraints on derived quantities like power dissipation, resonant frequency, or stiffness. Each constraint is a `CantileverMetricConstraint` referencing a `CantileverMetric` enum value.
+
+```python
+from piezod import (
+    CantileverMetric,
+    CantileverMetricConstraint,
+    optimize_performance,
+)
+
+constraints = [
+    CantileverMetricConstraint(CantileverMetric.POWER_DISSIPATION, maximum=2e-3),
+    CantileverMetricConstraint(CantileverMetric.OMEGA_VACUUM_HZ, minimum=5 * c.freq_max),
+    CantileverMetricConstraint(CantileverMetric.STIFFNESS, minimum=1e-3, maximum=1e1),
+]
+
+result = optimize_performance(
+    c,
+    force_resolution_goal(),
+    parameter_constraints={"max_v_bridge": 10.0},
+    metric_constraints=constraints,
+    n_starts=5,
+    max_iterations=10,
+    random_seed=0,
+)
+```
+
+`optimize_performance` runs SciPy's SLSQP (or L-BFGS-B when there are no nonlinear constraints) from random initial conditions and returns the best result once two starts agree within `convergence_tolerance` (1% by default), capped at `max_iterations` total runs. Use `optimize_performance_from_current` for a single-shot refinement of the existing design.
+
+```python
+print(f"force_resolution: {c.force_resolution() * 1e12:.1f} pN -> "
+      f"{result.optimized.force_resolution() * 1e12:.2f} pN")
+```
+
+Default geometric sanity constraints (`l/w >= 2`, `w/t >= 2`, `l_pr/w_pr >= 2`, `l_pr >= 2 um`) are added automatically; pass `default_aspect_constraints=False` to opt out.
+
+For a runnable end-to-end example with full output, see [`python/examples/optimization.py`](../python/examples/optimization.py). The example mirrors the MATLAB tutorial flow and typically improves Harley-1999-style force resolution by several hundred times.
+
+For implant-process-only optimization (geometry fixed, only `annealing_time/temp` and `implantation_energy/dose` varied), `CantileverImplantation.optimize_doping_for_hooge_noise` is also available with a Hooge-noise-limited default objective.
 
 ## Next Steps
 

python/examples/optimization.py

@@ -0,0 +1,93 @@
+"""Cantilever optimization example.
+
+Mirrors the optimization section of the MATLAB tutorial: an epitaxial
+piezoresistive cantilever is optimized for force resolution subject to
+power-dissipation and resonant-frequency constraints, with the bridge
+voltage capped at 10 V and the thickness floor at 1 um.
+
+Run:
+    uv run python examples/optimization.py
+"""
+
+from __future__ import annotations
+
+import warnings
+
+from piezod import (
+    CantileverEpitaxy,
+    CantileverMetric,
+    CantileverMetricConstraint,
+    force_resolution_goal,
+    optimize_performance,
+)
+
+
+def main() -> None:
+    # Starting design: 200 um x 20 um x 2 um epitaxial boron piezoresistor.
+    c = CantileverEpitaxy(
+        freq_min=1,
+        freq_max=1000,
+        l=200e-6,
+        w=20e-6,
+        t=2e-6,
+        l_pr_ratio=0.3,
+        v_bridge=2.0,
+        doping_type="boron",
+        dopant_concentration=1e19,
+        t_pr_ratio=0.3,
+    )
+    c.fluid = "vacuum"
+    c.number_of_piezoresistors = 4
+
+    print("=== Starting design ===")
+    print(f"  geometry: l={c.l * 1e6:.1f} um, w={c.w * 1e6:.1f} um, t={c.t * 1e6:.2f} um")
+    print(f"  v_bridge: {c.v_bridge:.2f} V")
+    print(f"  N_dopant: {c.dopant_concentration:.2e} cm^-3")
+    print(f"  t_pr_ratio: {c.t_pr_ratio:.3f}")
+    print(f"  force_resolution: {c.force_resolution() * 1e12:.2f} pN")
+    print(f"  power: {c.power_dissipation() * 1e3:.3f} mW")
+    print(f"  vacuum freq: {c.omega_vacuum_hz() / 1e3:.1f} kHz")
+    print(f"  stiffness: {c.stiffness():.3g} N/m")
+    print()
+
+    constraints = [
+        CantileverMetricConstraint(CantileverMetric.POWER_DISSIPATION, maximum=2e-3),
+        CantileverMetricConstraint(CantileverMetric.OMEGA_VACUUM_HZ, minimum=5 * 1000),
+        CantileverMetricConstraint(CantileverMetric.STIFFNESS, minimum=1e-3, maximum=1e1),
+    ]
+
+    print("Optimizing for force resolution with multi-start (this can take ~10 s)...")
+    with warnings.catch_warnings():
+        warnings.simplefilter("ignore")
+        result = optimize_performance(
+            c,
+            force_resolution_goal(),
+            parameter_constraints={"max_v_bridge": 10.0},
+            metric_constraints=constraints,
+            n_starts=5,
+            max_iterations=10,
+            random_seed=0,
+        )
+
+    opt = result.optimized
+    print(f"  ran {len(result.all_results)} starts; best objective {result.objective_value:.3f} pN")
+    print()
+
+    print("=== Optimized design ===")
+    print(f"  geometry: l={opt.l * 1e6:.1f} um, w={opt.w * 1e6:.1f} um, t={opt.t * 1e6:.3f} um")
+    print(f"  v_bridge: {opt.v_bridge:.2f} V")
+    print(f"  N_dopant: {opt.dopant_concentration:.2e} cm^-3")
+    print(f"  t_pr_ratio: {opt.t_pr_ratio:.3f}")
+    print(f"  force_resolution: {opt.force_resolution() * 1e12:.3f} pN")
+    print(f"  power: {opt.power_dissipation() * 1e3:.3f} mW (max {2.0:.1f})")
+    print(f"  vacuum freq: {opt.omega_vacuum_hz() / 1e3:.2f} kHz (min {5.0:.1f})")
+    print(f"  stiffness: {opt.stiffness():.3g} N/m")
+    print(f"  l/w ratio: {opt.l / opt.w:.2f} (default >= 2)")
+    print(f"  w/t ratio: {opt.w / opt.t:.2f} (default >= 2)")
+    print()
+    improvement = c.force_resolution() / opt.force_resolution()
+    print(f"Force resolution improved by {improvement:.1f}x")
+
+
+if __name__ == "__main__":
+    main()

python/src/piezod/__init__.py

@@ -20,6 +20,22 @@
     PiezoMaterial,
 )
 from piezod.cantilever_poly import CantileverPoly, Material
+from piezod.optimization import (
+    CantileverMetric,
+    CantileverMetricConstraint,
+    OptimizationResult,
+    StateVar,
+    charge_displacement_resolution_goal,
+    charge_force_resolution_goal,
+    displacement_resolution_goal,
+    force_noise_density_goal,
+    force_resolution_goal,
+    optimize_performance,
+    optimize_performance_from_current,
+    surface_stress_resolution_goal,
+    voltage_displacement_resolution_goal,
+    voltage_force_resolution_goal,
+)
 from piezod.piezoresistance import (
     CrystalOrientation,
     default_orientation,
@@ -33,6 +49,8 @@
     "CantileverDiffusion",
     "CantileverEpitaxy",
     "CantileverImplantation",
+    "CantileverMetric",
+    "CantileverMetricConstraint",
     "CantileverPiezoelectric",
     "CantileverPoly",
     "CrystalOrientation",
@@ -42,13 +60,25 @@
     "FluidType",
     "Material",
     "MetricConstraint",
+    "OptimizationResult",
     "PiezoMaterial",
     "PiezoresistorFromProfile",
+    "StateVar",
+    "charge_displacement_resolution_goal",
+    "charge_force_resolution_goal",
     "default_orientation",
     "diffusion_length_cm",
+    "displacement_resolution_goal",
+    "force_noise_density_goal",
+    "force_resolution_goal",
     "hooge_alpha_from_anneal",
+    "optimize_performance",
+    "optimize_performance_from_current",
     "pi_low_doping",
     "rotate_in_plane_stress",
+    "surface_stress_resolution_goal",
+    "voltage_displacement_resolution_goal",
+    "voltage_force_resolution_goal",
 ]
 
 __version__ = version("piezod")

python/src/piezod/cantilever.py

@@ -758,7 +758,9 @@ def print_performance(self) -> None:
         omega_damped_hz, Q = self.omega_damped_hz_and_Q()
         x, active_doping, total_doping = self.doping_profile()
         TMax_approx, TTip_approx = self.approxTempRise()
-        TMax, TTip = self.calculateMaxAndTipTemp()
+        # NOTE: the FD temperature solver (calculateMaxAndTipTemp /
+        # calculateTempProfile) has unported MATLAB-isms that raise on the
+        # default cantilever_type. Skip the F-D line until that is ported.
         thermoLimit = self.thermo_integrated() / self.force_sensitivity()
 
         print("=======================")
@@ -769,7 +771,9 @@ def print_performance(self) -> None:
         print(f"PR Length Ratio: {self.l_pr_ratio:g}")
 
         print(f"Force resolution (N): {self.force_resolution():g}")
-        print(f"Force noise at 1 kHz (fN): {self.force_noise_density(1e3):g}")
+        # force_noise_density returns a (1,1) array even for a scalar
+        # frequency input; squeeze to a Python float for the format spec.
+        print(f"Force noise at 1 kHz (fN): {float(np.squeeze(self.force_noise_density(1e3))) * 1e15:g}")
         print(f"Displacement resolution (m): {self.displacement_resolution():g}")
         print(f"Force sensitivity (V/N): {self.force_sensitivity():g}")
         print(f"Displacement sensitivity (V/m): {self.displacement_sensitivity():g}")
@@ -786,7 +790,6 @@ def print_performance(self) -> None:
         print(f"Sheet Resistance: {self.sheet_resistance()}")
         print(f"Power dissipation (mW): {self.power_dissipation() * 1e3:g}")
         print(f"Approx. Temp Rises (C) - Tip: {TTip_approx}  Max: {TMax_approx}")
-        print(f"F-D Temp Rises (C)     - Tip: {TTip}  Max: {TMax}")
 
         print(f"Integrated noise (V): {self.integrated_noise():g}")
         print(f"Integrated johnson noise (V): {self.johnson_integrated():g}")
@@ -802,7 +805,7 @@ def print_performance(self) -> None:
         print(f"Number of silicon resistors: {self.number_of_piezoresistors}")
         print(f"Si Thermal Conductivity (W/m-K): {self.k_base()}")
         print(f"E (GPa): {self.modulus() * 1e-9}")
-        print(f"Alpha: {self.alpha:g}")
+        print(f"Alpha: {self.alpha():g}")
 
         if self.cantilever_type == "step":
             print("=======================")
@@ -970,6 +973,7 @@ def hooge_integrated(self):
 
     # Johnson noise PSD from the entire Wheatstone bridge (V^2/Hz)
     def johnson_PSD(self, freq):
+        freq = np.atleast_1d(np.asarray(freq, dtype=float))
         resistance = self.resistance()  # resistance() includes contacts
         TPR = self.piezoresistor_temp()
         resistance *= 1 + TPR * self.TCR
@@ -1001,6 +1005,7 @@ def piezoresistor_temp(self):
     # Thermomechanical noise PSD
     # Units: V^2/Hz
     def thermo_PSD(self, freq):
+        freq = np.atleast_1d(np.asarray(freq, dtype=float))
         omega_damped_hz, Q_M = self.omega_damped_hz_and_Q()
         TPR = self.piezoresistor_temp()
         return (
@@ -1134,11 +1139,11 @@ def integrated_noise(self):
             + amplifier_integrated**2
         )
 
-    # Calculate the noise at a given frequency (V/rtHz)
+    # Calculate the noise at a given frequency (V/rtHz). Accepts either a
+    # scalar or an array; the result has the same shape semantics. np.sqrt
+    # is used (not math.sqrt) so array inputs work.
     def voltage_noise(self, freq):
-        return math.sqrt(
-            self.johnson_PSD(freq) + self.hooge_PSD(freq) + self.thermo_PSD(freq) + self.amplifier_PSD(freq)
-        )
+        return np.sqrt(self.johnson_PSD(freq) + self.hooge_PSD(freq) + self.thermo_PSD(freq) + self.amplifier_PSD(freq))
 
     def f_min_cumulative(self):
         frequency = np.logspace(math.log10(self.freq_min), math.log10(self.freq_max), Cantilever.numFrequencyPoints)
@@ -1747,10 +1752,12 @@ def calculateTempProfileTempDependent(self):
         T_final = T_current
         return x, Q, T_final
 
-    # Calculate the approximate thermal conductivity of the cantilever (W/m-K)
+    # Calculate the approximate thermal conductivity of the cantilever (W/m-K).
+    # k_x() already returns the scalar effective conductivity; the legacy
+    # ``k_x[0]`` indexing was a MATLAB-port leftover and raised on the float
+    # returned by the "approx" branch.
     def k_base(self):
-        k_x = self.k_x()
-        return k_x[0]
+        return self.k_x()
 
     # Calculate k(x) as a function of cantilever thickness
     # None/approx = Asheghi (1997) model
@@ -1884,9 +1891,9 @@ def force_noise_density(self, freq):
         return self.voltage_noise(freq) / self.force_sensitivity()
 
     def resonant_force_noise_density(self):
-        [omega_damped_hz, Q] = self.omega_damped_hz_and_Q()  ##ok<*NASGU>
-        freq = omega_damped_hz
-        return self.force_noise_density(freq) * 1e15
+        omega_damped_hz, _ = self.omega_damped_hz_and_Q()
+        psd = self.force_noise_density(np.atleast_1d(float(omega_damped_hz)))
+        return float(np.squeeze(psd)) * 1e15
 
     def surface_stress_resolution(self):
         return self.integrated_noise() / self.surface_stress_sensitivity()
@@ -2292,11 +2299,14 @@ def plot_noise_spectrum(self):
         freq = np.logspace(math.log10(self.freq_min), math.log10(self.freq_max), Cantilever.numFrequencyPoints)
 
         plt.figure()
-        plt.plot(freq, math.sqrt(self.johnson_PSD(freq)))
-        plt.plot(freq, math.sqrt(self.hooge_PSD(freq)))
-        plt.plot(freq, math.sqrt(self.thermo_PSD(freq)))
-        plt.plot(freq, math.sqrt(self.amplifier_PSD(freq)))
-        plt.plot(freq, self.voltage_noise(freq))
+        # np.sqrt (not math.sqrt) so each PSD array is rooted element-wise.
+        # Squeeze converts the (1, N) shape returned by johnson_PSD/thermo_PSD
+        # into (N,) so matplotlib gets matched 1-D x and y.
+        plt.plot(freq, np.squeeze(np.sqrt(self.johnson_PSD(freq))))
+        plt.plot(freq, np.squeeze(np.sqrt(self.hooge_PSD(freq))))
+        plt.plot(freq, np.squeeze(np.sqrt(self.thermo_PSD(freq))))
+        plt.plot(freq, np.squeeze(np.sqrt(self.amplifier_PSD(freq))))
+        plt.plot(freq, np.squeeze(self.voltage_noise(freq)))
         plt.gca().set_xscale("log")
         plt.gca().set_yscale("log")
         plt.ylabel("Noise Voltage Spectral Density (V/rtHz)")

python/src/piezod/cantilever_diffusion.py

@@ -17,6 +17,7 @@
 from scipy.special import erfc
 
 from piezod.cantilever import Cantilever
+from piezod.optimization.state import StateVar
 
 
 class CantileverDiffusion(Cantilever):
@@ -340,3 +341,21 @@ def doping_initial_conditions_random(self, parameter_constraints: dict[str, floa
         diffusion_temp_random = diffusion_temp_min + np.random.rand() * (diffusion_temp_max - diffusion_temp_min)
 
         return [diffusion_time_random, diffusion_temp_random]
+
+    def optimization_state_vars(self) -> tuple[StateVar, ...]:
+        """Declarative state spec for joint geometry + diffusion optimization.
+
+        Seven state variables: cantilever length, width, thickness,
+        piezoresistor length ratio, bridge bias voltage, diffusion time,
+        and diffusion temperature. Default bounds match MATLAB
+        ``cantileverDiffusion``.
+        """
+        return (
+            StateVar("l", 1e5, 10e-6, 3e-3),
+            StateVar("w", 1e7, 2e-6, 100e-6),
+            StateVar("t", 1e8, 1e-6, 100e-6),
+            StateVar("l_pr_ratio", 1e2, 0.01, 0.99),
+            StateVar("v_bridge", 1e1, 0.1, 10.0),
+            StateVar("diffusion_time", 1e-3, 5 * 60, 90 * 60),
+            StateVar("diffusion_temp", 1e-3, 273 + 800, 273 + 1000),
+        )