Empirical models of the thermospheric neutral density are routinely used by mission planners and systems engineers to perform orbit maintenance, collision avoidance, and estimate time and location of re-entry for spacecraft. These models have characteristic errors in neutral density below 10% during geomagnetic quiet time, but perform worse during intense geomagnetic activity, being unable to reproduce the significant increases in the neutral density that are observed during geomagnetic storms. Underestimation of the density during these conditions translates to errors in orbit propagation that reduce the accuracy of any resulting orbit predictions. These drawbacks directly translate into safety risks for astronauts and orbiting spacecraft, but also limit our understanding of the physics of neutral density enhancements. Numerous CubeSats with publicly available ephemeris in the form of two-line element (TLEs) sets orbit in this region. We present the Multifaceted Optimization Algorithm (MOA), a method to estimate the neutral density by minimizing the error between a modeled trajectory and a set of TLEs. Specifically, the algorithm estimates corrections to the inputs of the NRLMSISE-00 empirical density model, and applies those corrections along-track the SWARM spacecraft orbits. This results in orbit-averaged empirical densities below 10% error in magnitude, compared to errors in excess of 25\% for uncalibrated NLRMSISE-00.