Find non-negative roots of a system of equations using Scipy optimize.root

Question

I am trying to find the fixed points f(x)=x. I transform it to find the root of g(x)=f(x)-x.

from scipy.optimize import root
import numpy as np

def get_attraction(b1,b2, b3):
    c1 = 20 * b1 + 12 * b2 + 6 * b3
    c2 = 12 * b1 + 24 * b2 + 18 * b3
    c3 = 0 * b1 + 14 * b2 + 30 * b3
    return c1,c2,c3

def get_a_tau_max(p1_tau, p2_tau, p3_tau, a):
    b1_tau_l1_a1_a1 = p1_tau * (1 - a) + a
    b2_tau_l1_a1_a1 = p2_tau * (1 - a)
    b3_tau_l1_a1_a1 = p3_tau * (1 - a)

    a1_tau_l1_a1, a2_tau_l1_a1, a3_tau_l1_a1 = get_attraction(b1_tau_l1_a1_a1, b2_tau_l1_a1_a1,
                                                                   b3_tau_l1_a1_a1)


    b1_tau_l1_a1_a2 = p1_tau * (1 - a)
    b2_tau_l1_a1_a2 = p2_tau * (1 - a) + a
    b3_tau_l1_a1_a2 = p3_tau * (1 - a)
    a1_tau_l1_a2, a2_tau_l1_a2, a3_tau_l1_a2 = get_attraction(b1_tau_l1_a1_a2, b2_tau_l1_a1_a2,
                                                                   b3_tau_l1_a1_a2)

    b1_tau_l1_a1_a3 = p1_tau * (1 - a)
    b2_tau_l1_a1_a3 = p2_tau * (1 - a)
    b3_tau_l1_a1_a3 = p3_tau * (1 - a) + a
    a1_tau_l1_a3, a2_tau_l1_a3, a3_tau_l1_a3 = get_attraction(b1_tau_l1_a1_a3, b2_tau_l1_a1_a3,
                                                                   b3_tau_l1_a1_a3)

    a_tau_max = max(a1_tau_l1_a1, a2_tau_l1_a2, a3_tau_l1_a3)

    return a_tau_max

def get_prob(x, s1, s2, s3, a, lamda,p11,p12,p13,p21,p22,p23,p31,p32,p33):
    p1, p2, p3 = x[0], x[1], x[2]

    b1 = s1 * (1-a) + p1 * a
    b2 = s2 * (1-a) + p2 * a
    b3 = s3 * (1-a) + p3 * a

    c1,c2,c3=get_attraction(b1,b2,b3)
    
    t1 = get_a_tau_max(p11, p12, p13, a)
    t2 = get_a_tau_max(p21, p22, p23, a)
    t3 = get_a_tau_max(p31, p32, p33, a)

    a1 = c1+t1
    a2 = c2+t2
    a3 = c3+t3

    nom1 = np.exp(lamda*a1)
    nom2 = np.exp(lamda*a2)
    nom3 = np.exp(lamda*a3)
    dem = nom1 + nom2 + nom3

    p1_t = nom1 / dem
    p2_t = nom2 / dem
    p3_t = nom3 / dem

    return [p1_t - p1, p2_t - p2, p3_t - p3]


s1=2.0454331980638877e-42
s2=1.4352264582396677e-21
s3=1.0
a=11/13
lamda=4.0
p11=0.999999614695216
p12=3.853047896087925e-07
p13=5.242553502503789e-19
p21=1.425164081709055e-21
p22=0.9999999992759427
p23=7.240549367895844e-10
p31=2.03109266273484e-42
p32=1.4251640827409756e-21 
p33=1.0
result = root(get_prob, x0=np.array([0.1, 0.4, 0.5]), args=(s1, s2, s3, a, lamda,p11,p12,p13,p21,p22,p23,p31,p32,p33), method="lm")
print(result)

 message: The relative error between two consecutive iterates is at most 0.000000
 success: True
  status: 2
     fun: [ 1.773e-33  8.479e-28  0.000e+00]
       x: [-1.773e-33  7.462e-26  1.000e+00]
   cov_x: [[ 1.000e+00  0.000e+00  0.000e+00]
           [ 0.000e+00  1.000e+00 -3.309e-24]
           [ 0.000e+00 -3.309e-24  1.000e+00]]
    nfev: 13
    fjac: [[-1.000e+00 -3.309e-24 -1.000e+00]
           [-3.309e-24  1.000e+00 -3.309e-24]
           [ 0.000e+00  0.000e+00  1.000e+00]]
    ipvt: [3 2 1]
     qtf: [ 4.441e-16  3.888e-16  8.130e-22]

s1, s2, s3 are probabilities calculated using logitstic function, because it has nothing todo with a, I just supply the numbers in this example. Same for p11,p12,p13,and p21,p22,p23 and p31,p32,p33. They don't exactly sum up to 1.0, which I believe is because of precision (?).

I expect to get non-negative roots because of exp, depending on x0 I sometimes get a very small negative root.

Is this because of precision of the root method? How can I get non-negative roots?

I have seen some answers suggesting using optimize.minimize but this finding roots is already a part of my minimization problem.