Finding sets of binary sequences with low auto- and cross-correlation properties is a hard combinatorial optimization problem with numerous applications, including multiple-input-multiple-output (MIMO) radar and global navigation satellite systems (GNSS). The sum of squared correlations, also known as the integrated sidelobe level (ISL), is a quartic function in the variables and is a commonly-used metric of sequence set quality. In this paper, we show that the ISL minimization problem may be formulated as a mixed-integer quadratic program (MIQP). We then present a block coordinate descent (BCD) algorithm that iteratively optimizes over subsets of variables. The subset optimization subproblems are also MIQPs which may be handled more efficiently using specialized solvers than using brute-force search; this allows us to perform BCD over larger variable subsets than previously possible. Our approach was used to find sets of four binary sequences of lengths up to 1023 with better ISL performance than Gold codes and sequence sets found using existing BCD methods.