Calculating the yield limit Y c of a viscoplastic fluid flow is a challenging problem, often needing iteration in the rheological parameters to approach the plastic limit, as well as accurate computations that account properly for the yield stress and potentially adaptive meshing. For particle settling flows, in recent years calculating Y c has been accomplished analytically for many antiplane shear flow configurations and also computationally for many geometries, under either two dimensional (2D) or axisymmetric flow restrictions. Here we approach the problem of 3D particle settling and how to compute the yield limit directly, i.e. without iteratively changing the rheology to approach the plastic yield limit. The approach develops tools from optimization theory, taking advantage of the fact that Y c is defined via a minimization problem. We recast this minimization in terms of primal and dual variational problems, develop the necessary theory and finally implement a basic but workable algorithm. We benchmark results against accurate axisymmetric flow computations for cylinders and ellipsoids. This demonstrates the feasibility of directly computing Y c in multiple dimensions.