I am trying to maximize the solution of the system of algebraic equations for different values of parameters. The problem is related to the absorbing times in Markov chains.
Here is the code
Kk=10; (*chain length*)
ra = {1.5, 0.5}; \[Mu] = 0.5; \[Gamma]ab = {x, y}; \[Gamma]ba = {1.0 - x, 1.0 - y};
(* x and y here are variables that will be used in the maximization*)
Then I define the following lists, keeping p1,p2,x,y unspecified.
\[Omega]ip = Table[\[Mu] ra[[k]] i/(ra[[k]] i + (Kk - i)) (Kk - i)/Kk + (1 - \[Mu]) \[Gamma]ab[[k]] i/Kk (Kk - i)/Kk, {k, 1, 2}, {i, 1, Kk - 1}];
\[Omega]im = Table[\[Mu] (Kk - i)/(ra[[k]] i + (Kk - i)) i/Kk + (1 - \[Mu]) \[Gamma]ba[[k]] i/Kk (Kk - i)/Kk, {k, 1, 2}, {i, 1, Kk - 1}];
\[Gamma] = \[Omega]im/\[Omega]ip;
The system of algebraic equations are
eqns = Join[Table[v[i + 1] == \[Gamma][[1, i]] v[i] + 1/(1 - p2 - p1) 1/\[Omega]ip[[1, i]] p1 Sum[v[j] - vv[j], {j, 1, i}] - 1/\[Omega]ip[[1, i]], {i, 1, Kk - 1}],
Table[vv[i + 1] == \[Gamma][[2, i]] vv[i] + 1/(1 - p2 - p1) 1/\[Omega]ip[[2, i]] p2 Sum[vv[j] - v[j], {j, 1, i}] - 1/\[Omega]ip[[2, i]], {i, 1, Kk - 1}],
{Sum[v[i], {i, 1, Kk}] == 0}, {Sum[vv[i], {i, 1, Kk}] ==0}]
My goal is to maximize the unconditional absorbing time over x and y for given values of p1 and p2.
To do so I define a function
G[\[Pi]1_?NumericQ, \[Pi]2_?NumericQ] := Block[{ sol, v, Vv},
sol = Solve[eqns /. {p1 :> \[Pi]1, p2 :> \[Pi]2}, Join[Table[v[i], {i, 1, Kk}], Table[vv[i], {i, 1, Kk}]]];
v = First[v[1] /. sol]; Vv = First[vv[1] /. sol];
NMaximize[{1/(1 - \[Pi]2 - \[Pi]1) ((1 - \[Pi]2) v - \[Pi]1 Vv), 0 <= x && x <= 1 && 0 <= y && y <= 1}, {x, y}]] ;
However, it is too slow. For Kk=10 it takes 20 seconds for a pair of \[Pi]1 and [\Pi]2.
I need to do for Kk=100 and for 10000 pairs.
I tried to Solve eqns symbolically, and then insert the values of \[Pi]1 and \[Pi]2. It works for smaller Kk, but for larger Kk it does not.
Any help would be highly appreciated.
