Improving strategy to solve this constrained root-finding problem

Question

I am looking for roots of the function eqs defined below. It has $2k$ unknowns and returns a vector of $2k-1$ components.

The sole purpose of the following block is to define eqs, no need to try and understand it:

phi[p_, T_, t_] := -1/Sqrt[m[[p]]]*Sum[Cos[omegas[[j]]*t/2]/(omegas[[j]]*Sin[omegas[[j]]*T/2]) q[[p, j]]*q[[n,j]], {j, 1, n}];
psi[p_, sigma_, t_] :=  1/Sqrt[m[[p]]]* Sum[Sin[omegas[[j]]*t/2]/(Sin[omegas[[j]]*sigma/2]) q[[p, j]]* q[[n,j]], {j, 1, n}];
n=2;
m={.5,.5};
q={{-0.8506508083520401`,-0.5257311121191336`},{0.5257311121191336`,-0.8506508083520401`}};
omegas={3.23606797749979`,1.2360679774997898`};

aMat[sList_]:=Block[{tmp,T=Total[sList],k=Length[sList]},
tmp=Table[0,{i,1,k},{j,1,k}];
Table[tmp[[i,j]]=Sqrt[m[[n]]]*psi[n,T,2(Total[sList[[1;;j]]]-Total[sList[[1;;i]] ])-T],{i, 1, k}, {j, i+1,k}];
Table[tmp[[j,i]]=-tmp[[i,j]],{i,1,k},{j,i,k}];tmp
];

bMat[sList_]:=Block[{tmp,T=Total[sList],k=Length[sList]},
tmp=Table[0,{i,1,k},{j,1,k}];
Table[tmp[[i,j]]=phi[n,T,2(Total[sList[[1;;j]]]-Total[sList[[1;;i]] ])-T],{i, 1, k}, {j, i,k}];
Table[tmp[[j,i]]=tmp[[i,j]],{i,1,k},{j,i,k}];tmp
];

eqs[vars_]:=Block[{k=Length[vars]/2},(aMat[vars[[1;;k]]].vars[[k+1;;2k]])[[1;;-2]]~Join~(bMat[vars[[1;;k]]].vars[[k+1;;2k]]-ConstantArray[1,k])]

I am looking for solutions with all positive values and FindRoot cannot deal with constraints. Of course one choice could be to use FindRoot until the $2k-1$ components of solution are positive, but I think this is not an optimal strategy.

As suggested in this SE question, I wrote the problem as a constrained optimisation problem:

subdomain = Fold[And, Table[0 < vars[[i]], {i, 2, 2 k}]];
res = FindMinimum[{Total[eqs[vars0]^2], subdomain}, 
  Transpose[{vars[[2 ;; -1]], RandomReal[{1, 4}, 2 k - 1]}], 
  Method -> "InteriorPoint"]

I managed to get positive solutions for $k=6$ but it is not very robust and I would like to get solutions for larger $k$ (such as $15$, e.g.).

So my question is if someone see a way of improving my code or strategy in terms of robustness and efficiency.

A second optional question is to force the first $k$ values of the solution to be all different.

Note I have already asked a different question on the same set of equations here. Answers to the present question will probably make the previous one obsolete, so I'll delete it (or answer it if it is relevant).

Answer 1

One possibility is to add an equation that the sum of variables is some constant, say 1. You can actually enforce nonnegativity by augmenting the starting values, per refguide page for FindRoot. I won't repeat the lengthy setup code but just show the steps from there.

k = 9;
len = 2*k;
vars = Array[s, len];
eqns = Join[eqs[vars], {Total[vars] - 1}];
start = Thread[{vars, RandomReal[2, len]/len, 0, 1}];
FindRoot[eqns == 0, start]

During evaluation of In[59]:= FindRoot::reged: The point {0.03140964551874968,0.03874740609445434,0.0743937162451151,0.06308530530957773,0.1144968686002389,0.1358768793893159,0.00637481057165916,<<5>>,0.03178764934909824,0.1038665079060305,0.03872218273437743,0.01836388193562395,0.01507199254295023,0.} is at the edge of the search region {0.,1.} in coordinate 18 and the computed search direction points outside the region. >>

(* Out[59]= {s[1] -> 0.03140964551874968, s[2] -> 0.03874740609445434, 
 s[3] -> 0.0743937162451151, s[4] -> 0.06308530530957773, 
 s[5] -> 0.1144968686002389, s[6] -> 0.1358768793893159, 
 s[7] -> 0.00637481057165916, s[8] -> 0.01906482002753885, 
 s[9] -> 0.03327004567338578, s[10] -> 0.01749359129816518, 
 s[11] -> 0.0957320547892758, s[12] -> 0.05456770646648509, 
 s[13] -> 0.03178764934909824, s[14] -> 0.1038665079060305, 
 s[15] -> 0.03872218273437743, s[16] -> 0.01836388193562395, 
 s[17] -> 0.01507199254295023, s[18] -> 0.} *)

The message is not encouraging. My guess is it will be quite difficult to meet the constraints.