First, I would suggest that the strict inequalities in $0\leq x_1<x_2<\dots<x_N\leq1$ are changed into weak inequalities, such that we optimize over a closed set. (Then it is also clear that we can just optimize over the whole cube $[0,1]^N$: the vertices will we ordered in some way.)

Lets call the function we are trying to be optimized $f$. $$f=\sum_{S\in\binom{N}{K}} \prod_{i<j i,j \in S} (x_j - x_i)^2.$$

Maybe it is useful to note that we can assume $x_1=0$ and $x_N=1$, since otherwise we can just stretch the points to fit the unit interval, making $f$ larger.

I will give an answer for $N$ and $K$ small.

For $K=2$, the problem becomes a quadratic programming problem with $$f=x^TQx$$ for $$Q=\pmatrix{1&-1&-1&\dots&-1\\-1&1&-1&\dots &-1\\\dots&&&&\\-1&-1&-1&\dots&1}$$ Optimizing over the cube $[0,1]^N$ gives as solution, as mentioned in the comments, $$x_i\begin{cases}0 \text{ if }i\leq \left\lfloor\frac{N}{2}\right\rfloor\\1 \text{ otherwise }\end{cases}$$ The maximum attained is $$f=\left\lfloor\frac{N}{2}\right\rfloor\left\lceil\frac{N}{2}\right\rceil$$

Here are some pictures for $N=2$ and $N=3$.

The next image was made with sage by invoking the following command, which gives you also an animated 3d model

sum(implicit_plot3d((x-y)^2+(x-z)^2+(z-y)^2==.02+i*.3,(0,1),(0,1),(0,1),color=rainbow(4*7)[k-1-i]) for i in range(7)).show()

For $N=3$ the image shows level sets of $f$; they are cylinders arount the diagonal $\{(c,\dots,c) : c\in \mathbb{R}\}$. From that pictures it becomes clear that the optimizer must be the vertices of the cube that have the largest distance to the diagonal.

Here is a partially completed table of optimal solutions (value of $f$) for small $N$ and $K$ (Notice all values given are exact and not meant to be approximations): $$\begin{array}{lccccccc}K\backslash N&\bf2&\bf3&\bf4&\bf5&\bf6&\bf7&\bf8\\2&1&3&4&6&9&12&16\\3&-&\frac{1}{16}&\frac{1}{4}&&&\\4&-&-&\frac{1}{3125}&&&&\geq\frac{16}{3125}\end{array}$$

Let me list the optimal configuration of points. I omit $x_1=0$ and $x_N=1$. (Notice all values given are exact and not meant to be approximations)

$K=3$:

• $N=3$: $x_2=\frac{1}{2}$
• $N=4$: $x_2=x_3=\frac{1}{2}$
• $N=5$: $x_2=x_3=\frac{1}{2}$, $x_4=1$

$K=4$:

• $N=4$: $x_2=\frac{1}{2}-\frac{\sqrt{5}}{10}$, $x_3=\frac{1}{2}+\frac{\sqrt{5}}{10}$
• $N=8$ If we assume that the vertices come in pairs, as mentioned in the comments the points will be the same as in the case $N=4$. (Notice that the putatively optimal results mentioned in the comments $x_3=x_4=.275$ and $x_5=x_6=.725$ are almost but not quite optimal; better values would be $x_3=x_4=.276$ and $x_5=x_6=.724$ or even better $$x_3=x_4=0.27639320225002103035908263312687237645593816403885 \\x_5=x_6=0.72360679774997896964091736687312762354406183596115 $$

In all of these cases the points where values of $x_i$ lie are independent of $N$. Hence with the fact that they are the Fekete points for $N=K$, one could conjecture that they will also be the the Fekete points for other $N$ and search an optimum only among those.

First, I would suggest that the strict inequalities in $0\leq x_1<x_2<\dots<x_N\leq1$ are changed into weak inequalities:$0\leq x_1\leq x_2\leq\dots\leq x_N\leq1$, such that we optimize over a closed set. (Then it is also clear that we can just optimize over the whole cube $[0,1]^N$: the vertices will we ordered in some way.)

Lets call the function we are trying to be optimized $f$. $$f=\sum_{S\in\binom{N}{K}} \prod_{i<j i,j \in S} (x_j - x_i)^2.$$

Maybe it is useful to note that we can assume $x_1=0$ and $x_N=1$, since otherwise we can just stretch the points to fit the unit interval, making $f$ larger.

I will give an answer for $N$ and $K$ small.

For $K=2$, the problem becomes a quadratic programming problem with $$f=x^TQx$$ for $$Q=\pmatrix{N-1&-1&-1&\dots&-1\\-1&N-1&-1&\dots &-1\\\dots&&&&\\-1&-1&-1&\dots&N-1}$$ Optimizing over the cube $[0,1]^N$ gives as solution, as mentioned in the comments, $$x_i\begin{cases}0 \text{ if }i\leq \left\lfloor\frac{N}{2}\right\rfloor\\1 \text{ otherwise }\end{cases}$$ The maximum attained is $$f=\left\lfloor\frac{N}{2}\right\rfloor\left\lceil\frac{N}{2}\right\rceil$$

Here are some pictures for $N=2$ and $N=3$.

The next image was made with sage by invoking the following command, which gives you also an animated 3d model

sum(implicit_plot3d((x-y)^2+(x-z)^2+(z-y)^2==.02+i*.3,(0,1),(0,1),(0,1),color=rainbow(4*7)[k-1-i]) for i in range(7)).show()

For $N=3$ the image shows level sets of $f$; they are cylinders arount the diagonal $\{(c,\dots,c) : c\in \mathbb{R}\}$. From that pictures it becomes clear that the optimizer must be the vertices of the cube that have the largest distance to the diagonal.

Here is a partially completed table of optimal solutions (value of $f$) for small $N$ and $K$ obtained with the help of global optimization tools:

$$\begin{array}{lccccccc}K\backslash N&\bf2&\bf3&\bf4&\bf5&\bf6&\bf7&\bf8\\2&1&2&4&6&9&12&16\\3&-&\frac{1}{16}&\frac{1}{4}&\frac{1}{4}&\frac{1}{2}&\\4&-&-&\frac{1}{3125}&&&&\geq\frac{16}{3125}\end{array}$$

Let me list the optimal configuration of points. I omit $x_1=0$ and $x_N=1$. (Notice all values given are exact and not meant to be approximations)

$K=3$:

• $N=3$: $x_2=\frac{1}{2}$
• $N=4$: $x_2=x_3=\frac{1}{2}$
• $N=5$: $x_2=x_3=\frac{1}{2}$, $x_4=1$
• $N=6$: $x_2=0$, $x_3=x_4=\frac{1}{2}$, $x_5=1$

$K=4$:

• $N=4$: $x_2=\frac{1}{2}-\frac{\sqrt{5}}{10}$, $x_3=\frac{1}{2}+\frac{\sqrt{5}}{10}$
• $N=8$ If we assume that the vertices come in pairs, as mentioned in the comments the points will be the same as in the case $N=4$. (Notice that the putatively optimal results mentioned in the comments $x_3=x_4=.275$ and $x_5=x_6=.725$ are almost but not quite optimal; better values would be $x_3=x_4=.276$ and $x_5=x_6=.724$ or even better $$x_3=x_4=0.27639320225002103035908263312687237645593816403885 \\x_5=x_6=0.72360679774997896964091736687312762354406183596115 $$

In all of these cases the points where values of $x_i$ lie are independent of $N$. Hence with the fact that they are the Fekete points for $N=K$, one could conjecture that they will also be the the Fekete points for other $N$ and search an optimum only among those.