A new global optimization method is used to determine the distribution of earthquakes on a complex, connected fault system. The method, integer programming, has been advanced in the field of operations research, but has not been widely applied to geophysical problems until recently. In this application, we determine the optimal distribution of earthquakes on mapped faults to minimize the global misfit in slip rates for multi-fault ruptures. Integer programming solves for a decision vector composed of every possible location that a sample of earthquakes can occur on every fault, subject to slip-rate uncertainty constraints. Step over connections are straightforward to include, whereas branching fault connections are not. To include branching ruptures, we distinguish between individual multi-fault rupture paths, as opposed to formulating the integer-programming problem based on individual faults as in previous studies. The new method is applied to the complex fault system in the San Francisco Bay Area as a case study. Results from the integer-programming method are compared to those from a local optimization method, termed the greedy-sequential method. Several experiments using these two methods indicate that shape of the on-fault magnitude distributions and which branching faults are involved in multi-fault ruptures depend on how much emphasis is placed on fitting the target slip rate. In cases where the underlying data are not strong enough to warrant chasing the target slip rate, it is better to focus on the distribution of feasible results that better represents the uncertainty in the solutions imposed by the data.