trying to use Solve function to solve 4 linear equations and 4 variables (all probaibilites) but one of the resul returned is negative ? Am I using solve fn correctly

% In am interested in obtaining Steady state solution of 4 state ergodic
% continuous time markov chain where the long term probabilies in each
% state is given by Vector P = [ p1 p2 p3 p4] and 
% State Transition rate matrix A wehre 

% A=[-(2*lambda_System +lambda_com)          2*lambda_System       lambda_com  0;
%         mu       -(mu+lambda_standby+lambda_alert1) lambda_standby  lambda_alert1;
%          0            mu                   -(mu+lamda_alert2) lambda_alert2 0 0;
%          mu_repair    0                      0               -mu_repair;

% where lambda's and mu's are failure rates/ Recovery rates which are
% assumed as below based on Practical Field Data and are per hour 

% lambda_System = 0.00046 (Failure rate of System A assumed as once in 3 months)
% lambda_com = 0.00012 (Common mode failure rate as once in a year)
% lambda_alert1 = 3600 (alert1 to Station master rate for System A failure = 1 sec ie 3600/hr
% lambda_alert2 = 3600 (alter2 to Station master rate rate for Shut down = 1 sec ie 3600/hrto 
% mu  = 360 ( Auto reset recovery rate (once in 10 sec ie 360/hr)
% mu_repair 1 = 2.2 ( Minor repair recovery rate ( once in 30 minutes ie 2/hr) +
%                   (Shutdown or Major recovery rate (once in 5 hours ie 0.2/hr)
    
%     A =[-0.00104 0.00092 0.00012 0 ;
%     360       -360.00082 0.00046  3600;
%      0           0     -3600 3600;
%      2.2         0       0    - 2.2] 
%      

%   We obtain p1, p2 ...p4 by solving Equation 1 & 2 below:
%   P*A = 0   .... Equation 1 (Vectoe P and Matrix A are given above)
%   p1+p2+p3+4 = 1 ..... Equation 2
%  
%   we can obtain the values of p1 to p4 using the function solve

format short
syms eqn1 eqn2 eqn3  eqn4 eqn5 eqn6 eqn7 eqns sol p1sol p2sol p3sol p4sol p5sol p6sol   p1 p2 p3 p4 p5 p6
eqn1 = -p1*0.00104+p2*3600+ p4*2.2 == 0; 
eqn2 =  p1*0.00092 - p2*360.00082 == 0;  
eqn3 = p1*0.00012  + 0.00046*p2 - p3*3600 == 0; 
% eqn4 = p2*3600 + p3*3600 - p4*2.2 == 0;

% eqn5 = p1+p2+p3+p4==1; p4 = 1 - (p1+p2+p3)
% by substituting p4 from eqn5 in eqn4 we get eqn4 below

eqn4 = 2.2*p1 +3602.2*p2 + 3602.2*p3 == 2.2;

eqns = [eqn1, eqn2, eqn3, eqn4]
sol = solve([eqns],[p1,p2,p3,p4]);
p1sol = sol.p1
p2sol = sol.p2
p3sol = sol.p3
p4sol = sol.p4

double(sol.p1)
double(sol.p2)
double(sol.p3)
double(sol.p4)

% In am interested in obtaining Steady state solution of 4 state ergodic
% continuous time markov chain where the long term probabilies in each
% state is given by Vector P = [ p1 p2 p3 p4] and 
% State Transition rate matrix A wehre 

% A=[-(2*lambda_System +lambda_com)          2*lambda_System       lambda_com  0;
%         mu       -(mu+lambda_standby+lambda_alert1) lambda_standby  lambda_alert1;
%          0            mu                   -(mu+lamda_alert2) lambda_alert2 0 0;
%          mu_repair    0                      0               -mu_repair;

% where lambda's and mu's are failure rates/ Recovery rates which are
% assumed as below based on Practical Field Data and are per hour 

% lambda_System = 0.00046 (Failure rate of System A assumed as once in 3 months)
% lambda_com = 0.00012 (Common mode failure rate as once in a year)
% lambda_alert1 = 3600 (alert1 to Station master rate for System A failure = 1 sec ie 3600/hr
% lambda_alert2 = 3600 (alter2 to Station master rate rate for Shut down = 1 sec ie 3600/hrto 
% mu  = 360 ( Auto reset recovery rate (once in 10 sec ie 360/hr)
% mu_repair 1 = 2.2 ( Minor repair recovery rate ( once in 30 minutes ie 2/hr) +
%                   (Shutdown or Major recovery rate (once in 5 hours ie 0.2/hr)
    
%     A =[-0.00104 0.00092 0.00012 0 ;
%     360       -360.00082 0.00046  3600;
%      0           0     -3600 3600;
%      2.2         0       0    - 2.2] 
%      

%   We obtain p1, p2 ...p4 by solving Equation 1 & 2 below:
%   P*A = 0   .... Equation 1 (Vectoe P and Matrix A are given above)
%   p1+p2+p3+4 = 1 ..... Equation 2
%  
%   we can obtain the values of p1 to p4 using the function solve

format short
syms eqn1 eqn2 eqn3  eqn4 eqn5 eqn6 eqn7 eqns sol p1sol p2sol p3sol p4sol p5sol p6sol   p1 p2 p3 p4 p5 p6
eqn1 = -p1*0.00104+p2*3600+ p4*2.2 == 0; 
eqn2 =  p1*0.00092 - p2*360.00082 == 0;  
eqn3 = p1*0.00012  + 0.00046*p2 - p3*3600 == 0; 
% eqn4 = p2*3600 + p3*3600 - p4*2.2 == 0;

% eqn5 = p1+p2+p3+p4==1; p4 = 1 - (p1+p2+p3)
% by substituting p4 from eqn5 in eqn4 we get eqn4 below

eqn4 = 2.2*p1 +3602.2*p2 + 3602.2*p3 == 2.2;

eqns = [eqn1, eqn2, eqn3, eqn4]
sol = solve([eqns],[p1,p2,p3,p4]);
p1sol = sol.p1
p2sol = sol.p2
p3sol = sol.p3
p4sol = sol.p4

double(sol.p1)
double(sol.p2)
double(sol.p3)
double(sol.p4)


% In am interested in obtaining Steady state solution of 4 state ergodic
% continuous time markov chain where the long term probabilies in each
% state is given by Vector P = [ p1 p2 p3 p4] and 
% State Transition rate matrix A wehre 

% A=[-(2*lambda_System +lambda_com)          2*lambda_System       lambda_com  0;
%         mu       -(mu+lambda_standby+lambda_alert1) lambda_standby  lambda_alert1;
%          0            mu                   -(mu+lamda_alert2) lambda_alert2 0 0;
%          mu_repair    0                      0               -mu_repair;

% where lambda's and mu's are failure rates/ Recovery rates which are
% assumed as below based on Practical Field Data and are per hour 

% lambda_System = 0.00046 (Failure rate of System A assumed as once in 3 months)
% lambda_com = 0.00012 (Common mode failure rate as once in a year)
% lambda_alert1 = 3600 (alert1 to Station master rate for System A failure = 1 sec ie 3600/hr
% lambda_alert2 = 3600 (alter2 to Station master rate rate for Shut down = 1 sec ie 3600/hrto 
% mu  = 360 ( Auto reset recovery rate (once in 10 sec ie 360/hr)
% mu_repair 1 = 2.2 ( Minor repair recovery rate ( once in 30 minutes ie 2/hr) +
%                   (Shutdown or Major recovery rate (once in 5 hours ie 0.2/hr)
    
%     A =[-0.00104 0.00092 0.00012 0 ;
%     360       -360.00082 0.00046  3600;
%      0           0     -3600 3600;
%      2.2         0       0    - 2.2] 
%      

%   We obtain p1, p2 ...p4 by solving Equation 1 & 2 below:
%   P*A = 0   .... Equation 1 (Vectoe P and Matrix A are given above)
%   p1+p2+p3+4 = 1 ..... Equation 2
%  
%   we can obtain the values of p1 to p4 using the function solve

format short
syms eqn1 eqn2 eqn3  eqn4 eqn5 eqn6 eqn7 eqns sol p1sol p2sol p3sol p4sol p5sol p6sol   p1 p2 p3 p4 p5 p6
eqn1 = -p1*0.00104+p2*3600+ p4*2.2 == 0; 
eqn2 =  p1*0.00092 - p2*360.00082 == 0;  
eqn3 = p1*0.00012  + 0.00046*p2 - p3*3600 == 0; 
% eqn4 = p2*3600 + p3*3600 - p4*2.2 == 0;

% eqn5 = p1+p2+p3+p4==1; p4 = 1 - (p1+p2+p3)
% by substituting p4 from eqn5 in eqn4 we get eqn4 below

eqn4 = 2.2*p1 +3602.2*p2 + 3602.2*p3 == 2.2;

eqns = [eqn1, eqn2, eqn3, eqn4]
sol = solve([eqns],[p1,p2,p3,p4]);
p1sol = sol.p1
p2sol = sol.p2
p3sol = sol.p3
p4sol = sol.p4

double(sol.p1)
double(sol.p2)
double(sol.p3)
double(sol.p4)

Solve_ergodic_ctmc

eqns =

[3600*p2 - (13*p1)/12500 + (11*p4)/5 == 0, (23*p1)/25000 - (1583300350395579*p2)/4398046511104 == 0, (3*p1)/25000 + (23*p2)/50000 - 3600*p3 == 0, (11*p1)/5 + (18011*p2)/5 + (18011*p3)/5 == 11/5]

p1sol =

212971106914622718750000000/213873879410114629303113457

p2sol =

544258255749120000000/213873879410114629303113457

p3sol =

7099106441264547457/213873879410114629303113457

p4sol =

-789927167957101987500000/213873879410114629303113457

ans =

0.9958

ans =

2.5448e-06

ans =

3.3193e-08

ans =

-0.0037

>> c>>

Answer 1

[K Karthika]

0.9958 2.54486e-06 -0.00337 3.3193e-08

[K Karthika]

Usefull section

[K Karthika]

0.9958

[K Karthika]

2.5448e-06

[K Karthika]

3.3193e-08

Answer 2

[K Karthika]

-0.00337

Answer 3

[Walter Roberson]

% In am interested in obtaining Steady state solution of 4 state ergodic
% continuous time markov chain where the long term probabilies in each
% state is given by Vector P = [ p1 p2 p3 p4] and 
% State Transition rate matrix A wehre 

% A=[-(2*lambda_System +lambda_com)          2*lambda_System       lambda_com  0;
%         mu       -(mu+lambda_standby+lambda_alert1) lambda_standby  lambda_alert1;
%          0            mu                   -(mu+lamda_alert2) lambda_alert2 0 0;
%          mu_repair    0                      0               -mu_repair;

% where lambda's and mu's are failure rates/ Recovery rates which are
% assumed as below based on Practical Field Data and are per hour 

% lambda_System = 0.00046 (Failure rate of System A assumed as once in 3 months)
% lambda_com = 0.00012 (Common mode failure rate as once in a year)
% lambda_alert1 = 3600 (alert1 to Station master rate for System A failure = 1 sec ie 3600/hr
% lambda_alert2 = 3600 (alter2 to Station master rate rate for Shut down = 1 sec ie 3600/hrto 
% mu  = 360 ( Auto reset recovery rate (once in 10 sec ie 360/hr)
% mu_repair 1 = 2.2 ( Minor repair recovery rate ( once in 30 minutes ie 2/hr) +
%                   (Shutdown or Major recovery rate (once in 5 hours ie 0.2/hr)
    
%     A =[-0.00104 0.00092 0.00012 0 ;
%     360       -360.00082 0.00046  3600;
%      0           0     -3600 3600;
%      2.2         0       0    - 2.2] 
%      

%   We obtain p1, p2 ...p4 by solving Equation 1 & 2 below:
%   P*A = 0   .... Equation 1 (Vectoe P and Matrix A are given above)
%   p1+p2+p3+4 = 1 ..... Equation 2
%  
%   we can obtain the values of p1 to p4 using the function solve

format short
syms eqn1 eqn2 eqn3  eqn4 eqn5 eqn6 eqn7 eqns sol p1sol p2sol p3sol p4sol p5sol p6sol   p1 p2 p3 p4 p5 p6
eqn1 = -p1*0.00104+p2*3600+ p4*2.2 == 0; 
eqn2 =  p1*0.00092 - p2*360.00082 == 0;  
eqn3 = p1*0.00012  + 0.00046*p2 - p3*3600 == 0; 
% eqn4 = p2*3600 + p3*3600 - p4*2.2 == 0;

% eqn5 = p1+p2+p3+p4==1; p4 = 1 - (p1+p2+p3)
% by substituting p4 from eqn5 in eqn4 we get eqn4 below

eqn4 = 2.2*p1 +3602.2*p2 + 3602.2*p3 == 2.2;

eqns = [eqn1, eqn2, eqn3, eqn4]

eqns = 

sol1a = solve(eqns(1), p1)

sol1a = 

simplify(subs(eqns(2), p1, sol1a))

ans = 

sol1b = solve(eqns(2), p1)

sol1b = 

simplify(subs(eqns(1), p1, sol1b))

ans = 

Examine your equations.

If you solve the first one for p1 and substitute that into the second, then you end up with the sum of a positive value times p2 and a positive value times p4 equalling 0. That can only happen if one of p2 or p4 is negative.

If you solve the second one for p1 and substitute that into the first, then you end up with the exactly the same situation.

Therefore your equations cannot be solved without at least one of the p* values being negative.

[K Karthika]

35107874131426755p2+241892558110

[Walter Roberson]

That looks like a truncated version of what I posted.

In itself, if that has to equal zero, then

35107874131426755p2+241892558110 == 0
35107874131426755p2 == -241892558110
p2 = -241892558110 / 35107874131426755

so p2 would have to be negative.

When you have

A * p2 + B * p4 == 0
A * p2 == -B * p4
p2 = (-B/A) * p4

and if A and B are both positive then -B/A is negative, and it follows that p2 and p4 must be opposite signs (unless they are both zero)