How to minimize two functions at the same time with common variables

Question

I am trying to find the roots of three equations, and currently I am using minimize from scipy.

from scipy.optimize import basinhopping
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_prob(x, s1, s2, s3, a, lamda, t1, t2, t3):
    p1, p2 = x[0], x[1]
    p3 = 1 - p1 - p2

    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)

    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)**2+(p2_t - p2)**2+(p3_t - p3)**2

s1=1.0
s2=0.0
s3=0.0
a=1.0
lamda=-20
t1=30.0
t2=30.0
t3=30.0
result = basinhopping(get_prob, x0=[0.0, 0.25], niter=50, seed=np.random.seed(0),interval=10,
                              minimizer_kwargs={'args': (s1,s2,s3,a,lamda,t1,t2,t3),
                                                'method': "SLSQP",
                                                'tol': 1.0e-5,
                                                'bounds': [(0.0, 1.0), (0.0, 1.0)],
                                                'constraints': {'type': 'ineq',
                                                                'fun': lambda x: 1 - x[0] - x[1],
                                                                'jac': lambda x: np.full_like(x, -1)}})

print(result)

                    message: ['requested number of basinhopping iterations completed successfully']
                    success: True
                        fun: 1.646140610988733
                          x: [ 4.722e-02  9.926e-02]
                        nit: 50
      minimization_failures: 36
                       nfev: 2401
                       njev: 353
 lowest_optimization_result: message: Optimization terminated successfully
                             success: True
                              status: 0
                                 fun: 1.646140610988733
                                   x: [ 4.722e-02  9.926e-02]
                                 nit: 23
                                 jac: [-3.613e+00 -1.509e+00]
                                nfev: 164
                                njev: 23

Essentially, I want to find p1, p2, p3 (probability which should sum up to 1) such that p1=p1_t, p2=p2_t, p3=p3_t. Sometimes I got pretty large result.fun which is not so good for the root-finding problem.

I have also tried to use root and fixed_point. They perform worse than minimization.

Is there a way to minimize two functions at the time to improve the precision but without increase niter? Increasing number of iteration takes a long time, and it is not ideal for me.