Gradient system performance is increasingly constrained by Peripheral Nerve Stimulation (PNS). Nonetheless, gradient coil windings are optimized using the boundary element stream function method (BEM-SF) incorporating only non-biological metrics. We introduce direct incorporation of PNS constraints into a BEM-SF optimization in addition to the usual constraints. We pre-compute a novel PNS “oracle” matrix of each nerve segment’s likelihood to be excited by a given stream-function basis. The constraint is linear in current, so stimulation likelihood is a simple sum over all stream function bases. This allows convex optimization with a PNS constraint and examination of tradeoffs with linearity, inductance, torque and efficiency.