I think your formula is missing a logarithm of the difference of the longitudinal position of the two cross sections. That may not matter if your cross sections are taken at uniform distances from each other (but beware if they are not!). With the logarithm, your summation formula will be a first-order approximation to the integral found e.g. in these Stanford
Aircraft Design
course notes:
$$
F_{D}
\propto
\int
S^{\prime\prime}(x_1)
\,
S^{\prime\prime}(x_2)
\,
\ln \left| x_1 - x_2 \right|
\, \
dx_1
dx_2
$$
To minimize the drag, you have to vary the $S^{\prime\prime}$, naturally by varying the cross-sections $S$. And to do that, you need two things, a relation between $S$ and $S^{\prime\prime}$ (e.g. by choosing an analytic, parametric formula for the shape of S, and taking its second derivative) and a constraint that prevents it from going to zero (e.g. a constant volume or a constant maximum cross section).
If you go with my suggestion above (analytic cross section formula and the most straight-forward cross-sectional constraint), you may find that you can solve the problem analytically (finding a Sears-Haack body). Of course an approximate numerical solution is always possible, and may be convenient or necessary depending on your constraints (and formula for linking the cross sections to their second derivative).
