Findroot apparently claiming to find singular Jacobians at all initial points for a system of 19 equations with 19 unknowns

Question

The following code symbolically computes a system of equations and then applies FindRoot to it. When I run code the message 'Encountered a singular Jacobian...' comes up, for all initial points that I have tried.

k = 0.5;
fA = Max[0, Min[1, (1/h)*(Log[(pAC - k)/(pA - k)] + pAC/(pAC - k) - pA*(pA - k))]]; 
fB = Max[0, Min[1, (1/h)*(Log[(pBC - k)/(pB - k)] + pBC/(pBC - k) - pB*(pB - k))]]; 
zAB = Max[0, Min[1, 1/2 + (1/(2*g))*(Log[(pB - k)/(pA - k)] + pB/(pB - k) - pA/(pA - k))]]; 
zC = Max[0, Min[1, 1/2 + (1/(2*g))*(Log[(pBC - k)/(pAC - k)] + pBC/(pBC - k) - pAC/(pAC - k))]]; 
rAA = fA*zAB - (1/2)*Max[0, zAB - zC]*Max[0, fA - fB]; 
rBB = fB*(1 - zAB) - (1/2)*Max[0, (1 - zAB) - (1 - zC)]*Max[0, fB - fA]; 
rAC = (1 - fA)*zC - (1/2)*Max[0, zC - zAB]*Max[0, (1 - fA) - (1 - fB)]; 
rBC = (1 - fB)*(1 - zC) - (1/2)*Max[0, (1 - zC) - (1 - zAB)]*Max[0, (1 - fB) - (1 - fA)]; 
disutilityAA = fA*zAB*(h*(fA/2) + g*(zAB/2)) - (1/2)*Max[0, zAB - zC]*Max[0, fA - fB]*(h*(fA - (1/3)*(fA - fB)) + g*(zAB - (1/3)*(zAB - zC))); 
disutilityBB = fB*(1 - zAB)*(h*(fB/2) + g*((1 - zAB)/2)) - (1/2)*Max[0, (1 - zAB) - (1 - zC)]*Max[0, fB - fA]*
     (h*(fB - (1/3)*(fB - fA)) + g*((1 - zAB) - (1/3)*((1 - zAB) - (1 - zC)))); 
disutilityAC = (1 - fA)*zC*(h*((1 - fA)/2) + g*(zC/2)) - (1/2)*Max[0, zC - zAB]*Max[0, (1 - fA) - (1 - fB)]*
     (h*((1 - fA) - (1/3)*((1 - fA) - (1 - fB)) + g*(zC - (1/3)*(zC - zAB)))); 
disutilityBC = (1 - fB)*(1 - zC)*(h*((1 - fB)/2) + g*((1 - zC)/2)) - (1/2)*Max[0, (1 - zC) - (1 - zAB)]*(h*((1 - fB) - (1/3)*((1 - fB) - (1 - fA))) + 
      g*((1 - zC) - (1/3)*((1 - zC) - (1 - zAB))));
everythingintermsofplayersactions = {xA -> 1/(pA - k), xAC -> 1/(pAC - k), xB -> 1/(pB - k), xBC -> 1/(pBC - k)};
\[Pi]A = (pA/(1 + tA) - c)*xA*(mA + rAA*mC) + (pAC/(1 + tC) - c)*xAC*rAC*mC //. everythingintermsofplayersactions; 
\[Pi]B = (pB/(1 + tB) - c)*xB*(mB + rBB*mC) + (pBC/(1 + tC) - c)*xBC*rBC*mC //. everythingintermsofplayersactions; 
focA = D[\[Pi]A, pA] == 0; 
focAC = D[\[Pi]A, pAC] == 0; 
focB = D[\[Pi]B, pB] == 0; 
focBC = D[\[Pi]B, pBC] == 0; 
dependenceontaxes = {pA -> pA[tA, tB, tC], pB -> pB[tA, tB, tC]}; 
focsp = {focA, focAC, focB, focBC} /. dependenceontaxes; 
t = {tA, tB, tC}; 
dfocspbydt = D[focsp, {t}]; 
\[Pi]Aasafunctionoft = \[Pi]A /. dependenceontaxes; 
\[Pi]Basafunctionoft = \[Pi]B /. dependenceontaxes; 
sA = \[Pi]A + (tA/(1 + tA))*pA*xA*(mA + mC*rAA) + mA*(xA*k - xA*pA) //. everythingintermsofplayersactions /. dependenceontaxes; 
sB = \[Pi]B + (tB/(1 + tB))*pB*xB*(mB + mC*rBB) + mB*(yB*k - xB*pB) //. everythingintermsofplayersactions /. dependenceontaxes; 
  sC = (tC/(1 + tC))*mC*(rBC*xBC*pBC + rAC*xAC*pAC) + mC*(k*(rAA*xA + rBB*xB + rAC*xAC + rBC*xBC) - rAA*xA*pA - rBB*xB*pB - rAC*xAC*pAC - rBC*xBC*pBC) - 
      disutilityAA - disutilityAC - disutilityBB - disutilityBC //. everythingintermsofplayersactions /. dependenceontaxes; 
focsA = D[sA, tA] == 0; 
focsB = D[sB, tB] == 0; 
focsC = D[sC, tC] == 0; 
allconditions = Join[Flatten[focsp], Flatten[dfocspbydt], {focsA}, {focsB}, {focsC}]; 
allconditionswithoutfunctions = 
   allconditions /. pAC[tA, tB, tC] -> pAC /. pA[tA, tB, tC] -> pA /. pB[tA, tB, tC] -> pB /. pBC[tA, tB, tC] -> pBC /. Derivative[1, 0, 0][pA][tA, tB, tC] -> 
                dpAdtA /. Derivative[0, 1, 0][pA][tA, tB, tC] -> dpAdtB /. Derivative[0, 0, 1][pA][tA, tB, tC] -> dpAdtC /. Derivative[1, 0, 0][pAC][tA, tB, tC] -> 
             dpACdtA /. Derivative[0, 1, 0][pAC][tA, tB, tC] -> dpACdtB /. Derivative[0, 0, 1][pAC][tA, tB, tC] -> dpACdtC /. 
         Derivative[1, 0, 0][pB][tA, tB, tC] -> dpBdtA /. Derivative[0, 1, 0][pB][tA, tB, tC] -> dpBdtB /. Derivative[0, 0, 1][pB][tA, tB, tC] -> dpBdtC /. 
      Derivative[1, 0, 0][pBC][tA, tB, tC] -> dpBCdtA /. Derivative[0, 1, 0][pBC][tA, tB, tC] -> dpBCdtB /. Derivative[0, 0, 1][pBC][tA, tB, tC] -> dpBCdtC; 
parameters = {\[Rho] -> 1/2, c -> 1, mA -> 1/3, mB -> 1/3, mC -> 1/3, g -> 0.15, h -> 0.2}; 
allconditionswithoutfunctionswithnumbers = allconditionswithoutfunctions /. parameters; 
FindRoot[{allconditionswithoutfunctionswithnumbers}, {{tA, 0.43}, {tB, 0.44}, {tC, 0.88}, {pAC, 6.1}, {pA, 4.}, {pBC, 7.1}, {pB, 3.2}, {dpAdtA, 1.1}, {dpAdtB, 1.2}, 
   {dpAdtC, 1.05}, {dpACdtA, 1.06}, {dpACdtB, 1.03}, {dpACdtC, 1.02}, {dpBdtA, 1.04}, {dpBdtB, 1.03}, {dpBdtC, 1.01}, {dpBCdtA, 1.02}, {dpBCdtB, 1.04}, 
   {dpBCdtC, 1.03}}]

The system consists of 19 equations and has 19 variables. Is it possible that this is due to a problem with FindRoot? (From the context of the problem it seems that there should be a solution and that the solutions should be locally unique.)

When I change the value of k to k=0 at the beginning, then I instead always get the error message 'The line search decreased the step size to within tolerance specified \ by AccuracyGoal and PrecisionGoal but was unable to find a sufficient \ decrease in the merit function. You may need more than \ MachinePrecision digits of working precision to meet these tolerances.'

Answer 1

In the meantime, I did some experiments. I repaced all occurences of Max[0,x] by Ramp[x] and Min[1,x] by 1 - Ramp[1 - x] since these can be easily differentiated by Mathematica.

Moreover, I computed the derivative of your expression

Off[Part::partd];
vars = {tA, tB, tC, pAC, pA, pBC, pB, dpAdtA, dpAdtB, dpAdtC, dpACdtA,
    dpACdtB, dpACdtC, dpBdtA, dpBdtB, dpBdtC, dpBCdtA, dpBCdtB, 
   dpBCdtC};
XX = Table[X[[i]], {i, 1, Length[vars]}];
expr = allconditionswithoutfunctionswithnumbers[[All, 1]];
F = X \[Function] Evaluate[N[expr /. Thread[vars -> XX]]];
DF = X \[Function] Evaluate[N[D[F[XX], {XX, 1}]]];

Let's calculate the derivative at the initial point and check its rank:

X0 = {0.43, 0.44, 0.88, 6.1, 4., 7.1, 3.2, 1.1, 1.2, 1.05, 1.06, 1.03,
    1.02, 1.04, 1.03, 1.01, 1.02, 1.04, 1.03};
Dimensions[DF[X0]]
MatrixRank[DF[X0]]

{19, 19}

12

Nope, it's not invertible, even with lots of good will.