The principle of finite descent

Preliminaries

Select 10 points at random in a unit square :

n=10;

red=RandomReal[{0,1},{n,2},WorkingPrecision->12];   rouge=red;
blue=RandomReal[{0,1},{n,2},WorkingPrecision->12];    bleu=blue;

Permute elements i and j in a list :

liste={{a1,a2},{b1,b2},{c1,c2},{d1,d2}}

exch[liste_,i_,j_]:=Flatten[{z=liste[[j]];liste/.{liste[[j]]->liste[[i]],liste[[i]]->z}},1]

exch[liste,3,2]

Test if segments i (joining red[[i]] and blue[[i]]) and j (joining red[[j]] and blue[[j]]) intersect :

test[{i_,j_}]:=((Det[{rouge3D[[i]],bleu3D[[j]],rouge3D[[j]]}] Det[{bleu3D[[i]],bleu3D[[j]],rouge3D[[j]]}])<0 && (Det[{rouge3D[[j]],bleu3D[[i]],rouge3D[[i]]}] Det[{bleu3D[[j]],bleu3D[[i]],rouge3D[[i]]}])<0)

(See Annexe below)

Finite descent (first example)

10 red points and 10 blues points in a unit square :

n=10;

A pleasant random choice constantly valid in this file :

red={{0.0053798,0.33624394},{0.423924225,0.056597131},{0.111222888,0.553108443},{0.2178671076,0.18560035671},{0.40593609037,0.04100032948},{0.58612420387,0.62245760504},{0.9425556202,0.1233684163},{0.199335730108,0.4698589344},{0.04701523255,0.42571309118},{0.87918160556,0.17446781699}};

blue={{0.26195419024,0.5609855375},{0.8041082577,0.488814573},{0.0673532619,0.210971476929},{0.30939673809,0.6669391066},{0.4888056944,0.9407516235},{0.568572576611,0.55565696683},{0.22113230546,0.93019714635},{0.82160565233061,0.5576455827140},{0.20638563594,0.1806507434},{0.295234104635,0.24453300739}};

Initial display :

The 20 points :

rouge=red;rouge3D=Append[#,1]&/@rouge ;    
bleu=blue;bleu3D=Append[#,1]&/@bleu;

Graphics[{PointSize[Medium],Line[{{0,0},{1,0},{1,1},{0,1},{0,0}}],Red,Point[rouge],Blue,Point[bleu]}]

One-to-one pairing of the red and blue points (In the order 1-1, 2-2, ..., 10-10) (Total length, L,  of the straight segments displayed) :

Graphics[{PointSize[Medium],Line[{{0,0},{1,0},{1,1},{0,1},{0,0}}],Table[Line[{rouge[[k]],bleu[[k]]}],{k,n}],Red,Point[rouge],Blue,Point[bleu],Black,Inset[Total[Table[ArcLength[Line[{rouge[[k]],bleu[[k]]}]],{k,n}]],{0.5,-0.04}]}]

Progressive uncrossing of the segments. If the final  uncrossing is incomplete, the cell must be manually reactivated until complete uncrossing is achieved :

Do[{If[test[{i,j}],{bleu=exch[bleu,i,j],bleu3D=exch[bleu3D,i,j],graphes=Append[graphes,Graphics[{PointSize[Medium],Line[{{0,0},{1,0},{1,1},{0,1},{0,0}}],Table[Line[{rouge[[k]],bleu[[k]]}],{k,n}],Red,Point[rouge],Blue,Point[bleu],Black,Inset[Total[Table[ArcLength[Line[{rouge[[k]],bleu[[k]]}]],{k,n}]],{0.5,-0.04}]}]]}, ]},{i,1,n},{j,i+1,n}];
graphes

GraphicsGrid[Partition[graphes,UpTo[4]]]

Finite descent (second example)

The previous calculation started with bleu=Permutations[blue][[1]] (=blue i.e. Identity permutation). One can do better if starting  with bleu=Permutations[blue][[1066281]] :

Position[Permutations[{1,2,3,4,5,6,7,8,9,10}],{3,10,5,6,9,7,2,4,1,8}]

Permutations[{1,2,3,4,5,6,7,8,9,10}][[1066281]]

n=10;k=1066281;

bleu=Permutations[blue][[k]]

bleu3D=Append[#,1]&/@bleu;

New starting point :

graphes={Graphics[{PointSize[Medium],Line[{{0,0},{1,0},{1,1},{0,1},{0,0}}],Table[Line[{rouge[[k]],bleu[[k]]}],{k,n}],Red,Point[rouge],Blue,Point[bleu],Black,Inset[Total[Table[ArcLength[Line[{rouge[[k]],bleu[[k]]}]],{k,n}]],{0.5,-0.04}]}]}

Do[{If[test[{i,j}],{bleu=exch[bleu,i,j],bleu3D=exch[bleu3D,i,j],graphes=Append[graphes,Graphics[{PointSize[Medium],Line[{{0,0},{1,0},{1,1},{0,1},{0,0}}],Table[Line[{rouge[[k]],bleu[[k]]}],{k,n}],Red,Point[rouge],Blue,Point[bleu],Black,Inset[Total[Table[ArcLength[Line[{rouge[[k]],bleu[[k]]}]],{k,n}]],{0.5,-0.04}]}]]}, ]},{i,1,n},{j,i+1,n}];
graphes

GraphicsGrid[Partition[graphes,UpTo[4]]]

Annexe : Criterion that decides the intersection between two straight segments, AB and CD

Intersection point (X = {xx, yx}) between two straight lines, ab et cd :

resu=Flatten[Solve[{((x-xa)/(xb-xa))==((y-ya)/(yb-ya)),((x-xc)/(xd-xc))==((y-yc)/(yd-yc))},{x,y}]]//FullSimplify

{xx,yx}={x,y}/.resu

{xd-xx,yd-yx}.{xc-xx,yc-yx}//Simplify

residudc=(xc (-ya+yb)+xb (ya-yc)+xa (-yb+yc))(xd (-ya+yb)+xb (ya-yd)+xa (-yb+yd))      (*Must be negative if segments ab and cd crosses *)

{xb-xx,yb-yx}.{xa-xx,ya-yx}//Simplify

residuba=(xd (ya-yc)+xa (yc-yd)+xc (-ya+yd)) (xd (yb-yc)+xb (yc-yd)+xc (-yb+yd))       (*Must be negative if segments ab and cd crosses *)

Det[{{xb,yb,1},{xc,yc,1},{xd,yd,1}}] Det[{{xa,ya,1},{xc,yc,1},{xd,yd,1}}]==residuba //Simplify

Det[{{xc,yc,1},{xb,yb,1},{xa,ya,1}}] Det[{{xd,yd,1},{xb,yb,1},{xa,ya,1}}]==residudc //Simplify

((Det[{{xb,yb,1},{xc,yc,1},{xd,yd,1}}//MatrixForm]Det[{{xa,ya,1},{xc,yc,1},{xd,yd,1}}//MatrixForm]<0)/.{xa->(Subscript[x, A]),xb->(Subscript[x, B]),xc->(Subscript[x, C]),xd->(Subscript[x, D]),ya->(Subscript[y, A]),yb->(Subscript[y, B]),yc->(Subscript[y, C]),yd->(Subscript[y, D])}) && ((Det[{{xc,yc,1},{xb,yb,1},{xa,ya,1}}//MatrixForm]Det[{{xd,yd,1},{xb,yb,1},{xa,ya,1}}//MatrixForm]<0)/.{xa->(Subscript[x, A]),xb->(Subscript[x, B]),xc->(Subscript[x, C]),xd->(Subscript[x, D]),ya->(Subscript[y, A]),yb->(Subscript[y, B]),yc->(Subscript[y, C]),yd->(Subscript[y, D])})