AbstractEmpirical models of the thermospheric density are routinely used to perform orbit maintenance, satellite collision avoidance, and estimate time and location of re‐entry for spacecraft. These models have characteristic errors in the thermospheric density below 10% during geomagnetic quiet time but are unable to reproduce the significant increase and subsequent recovery in the density 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 risk the safety of astronauts and orbiting spacecraft and also limit understanding of the physics of thermospheric 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 thermospheric density by minimizing the error between a modeled trajectory and a set of TLEs. The algorithm first estimates a representative cross‐sectional area for several reference CubeSats during the quiet time 3 weeks prior to the storm, and then estimates modifications to the inputs of the NRLMSISE‐00 empirical density model in order to minimize the difference between the modeled and TLE‐provided semimajor axis of the CubeSats. For validation, the median value of the modifications across all CubeSats are applied along the Swarm spacecraft orbits. This results in orbit‐averaged empirical densities below 10% error in magnitude during a geomagnetic storm, compared to errors in excess of 25% for uncalibrated NLRMSISE‐00 when compared to Swarm GPS‐derived densities.