Can't NSolve system of exponential and linear equations in several variables

Question

At the moment I am stuck with the following code, when specified to use root , quadratic or log utility functions all it fine, only exponentials like:

$\small \left\{-e^{-c_{0,1}}-0.25 e^{-c_{1,1}}-0.25 e^{-c_{2,1}}-0.5 e^{-c_{3,1}}\right\}$ s.t $ \small\left\{\text{pr}_1 \left(c_{1,1}-5\right)+\text{pr}_2 \left(c_{2,1}-10\right)+\text{pr}_3 \left(c_{3,1}-10\right)+c_{0,1}-10\right\}$

won't give out any solution, even when all equations (market-clearing like sum of agents consumption equal to endowments of state) are specified.

The whole file can be downloaded here. Every section is explained. The utility-function is modified from-to:

    Subscript[Vf, 1][j_, i_, \[Alpha]_] := consall[[i, j]] - \[Alpha]*consall[[i, j]]^2
    Subscript[Vf, 1][j_, i_, \[Alpha]_] := -Exp[-consall[[i, j]]]

First one gives output of the two solutions:

$$\small \left\{\left\{c_{-1,1}\to 0.816667,c_{-1,2}\to 0.883333,c_{0,1}\to 9.16667,c_{0,2}\to 5.83333,c_{1,1}\to 9.16667,c_{1,2}\to 5.83333,c_{2,1}\to 9.16667,c_{2,2}\to 5.83333,c_{3,1}\to 50.,c_{3,2}\to 50.,\text{pr}_1\to 0.25,\text{pr}_2\to 0.25,\text{pr}_3\to 0.\right\},\left\{c_{-1,1}\to 0.824265,c_{-1,2}\to 0.875735,c_{0,1}\to 8.78673,c_{0,2}\to 6.21327,c_{1,1}\to 8.78673,c_{1,2}\to 6.21327,c_{2,1}\to 8.78673,c_{2,2}\to 6.21327,c_{3,1}\to 11.211,c_{3,2}\to 8.78896,\text{pr}_1\to 0.25,\text{pr}_2\to 0.25,\text{pr}_3\to 0.470588\right\}\right\}$$

Exponential only yield the message: "Solve was unable to solve the system with inexact coefficients or the system obtained by direct rationalization of inexact numbers present in the system. Since many of the methods used by Solve require exact input, providing Solve with an exact version of the system may help."

All parameters are given numerical(random or specified) values so that it should solve the system. FindRoot is not suitable as the number of solutions are too big to find the right one even if ranges are specified (ranges of mx(n-1)+n-1 variables would have to be stated !!) Any help would be great, Thanks!

Answer 1

I'm not going to give a full answer because it will appear as a mess of Subscript notation (if you want to avoid such things, use a different notation at least until you have a debugged code base). There are at least a couple of routes one might take to get solutions. One, which I did with success, is to take the solutions to the original system as inputs to FindRoot for the revised system. From this I got two distinct solutions.

tester = vartab /. Solve[Flatten[fefe2], vartab];

eqns3

(* Out[346]= {E^-Subscript[c, 0, 1] == Subscript[c, -1, 1], 
 0.25 E^-Subscript[c, 1, 1] == Subscript[pr, 1] Subscript[c, -1, 1], 
 0.25 E^-Subscript[c, 2, 1] == Subscript[pr, 2] Subscript[c, -1, 1], 
 0.5 E^-Subscript[c, 3, 1] == 
  Subscript[pr, 3] Subscript[c, -1, 1], -10 + Subscript[c, 0, 1] + 
   Subscript[pr, 1] (-5 + Subscript[c, 1, 1]) + 
   Subscript[pr, 2] (-10 + Subscript[c, 2, 1]) + 
   Subscript[pr, 3] (-10 + Subscript[c, 3, 1]) == 0, 
 E^-Subscript[c, 0, 2] == Subscript[c, -1, 2], 
 0.25 E^-Subscript[c, 1, 2] == Subscript[pr, 1] Subscript[c, -1, 2], 
 0.25 E^-Subscript[c, 2, 2] == Subscript[pr, 2] Subscript[c, -1, 2], 
 0.5 E^-Subscript[c, 3, 2] == 
  Subscript[pr, 3] Subscript[c, -1, 2], -5 + Subscript[c, 0, 2] + 
   Subscript[pr, 1] (-10 + Subscript[c, 1, 2]) + 
   Subscript[pr, 2] (-5 + Subscript[c, 2, 2]) + 
   Subscript[pr, 3] (-10 + Subscript[c, 3, 2]) == 
  0, -15 + Subscript[c, 0, 1] + Subscript[c, 0, 2] == 
  0, -15 + Subscript[c, 1, 1] + Subscript[c, 1, 2] == 
  0, -15 + Subscript[c, 2, 1] + Subscript[c, 2, 2] == 0} *)

start1 = Transpose[{vartab, tester[[1]]}]
rt1 = FindRoot[eqns3, start1, MaxIterations -> 1000]

(* Out[342]= {{Subscript[c, -1, 1], 0.816666666667}, {Subscript[c, -1, 
  2], 0.883333333333}, {Subscript[c, 0, 1], 
  9.16666666667}, {Subscript[c, 0, 2], 5.83333333333}, {Subscript[c, 
  1, 1], 9.16666666667}, {Subscript[c, 1, 2], 
  5.83333333333}, {Subscript[c, 2, 1], 9.16666666667}, {Subscript[c, 
  2, 2], 5.83333333333}, {Subscript[c, 3, 1], 50.}, {Subscript[c, 3, 
  2], 50.}, {Subscript[pr, 1], 0.25}, {Subscript[pr, 2], 
  0.25}, {Subscript[pr, 3], 0.}} *)

(* Out[343]= {Subscript[c, -1, 1] -> 0.000104464143832, 
 Subscript[c, -1, 2] -> 0.00292829969482, 
 Subscript[c, 0, 1] -> 9.16666666667, 
 Subscript[c, 0, 2] -> 5.83333333333, 
 Subscript[c, 1, 1] -> 9.16666666667, 
 Subscript[c, 1, 2] -> 5.83333333333, 
 Subscript[c, 2, 1] -> 9.16666666667, 
 Subscript[c, 2, 2] -> 5.83333333333, Subscript[c, 3, 1] -> 50., 
 Subscript[c, 3, 2] -> 50., Subscript[pr, 1] -> 0.25, 
 Subscript[pr, 2] -> 0.25, Subscript[pr, 3] -> 6.56496072263*10^-54} *)

Similar for the other one.

The second possibility is to separate into polynomial and nonpolynomial equations, solve for an appropriate set of variables in the polynomial subsystem, and substitute those into the remaining equations. One way to go about this is indicated below.

eqns4 = Select[eqns3, FreeQ[#, E] &];
solns = Solve[
  eqns4, {Subscript[c, 0, 1], Subscript[pr, 2], Subscript[c, 2, 2], 
   Subscript[c, 1, 2], Subscript[c, 3, 2]}]
eqns5 = Together[eqns3 /. solns] /. True -> Sequence[]
Solve[eqns5[[1]]]

It seems to be quite slow at the last step though and I do not know if it will complete in reasonable time.