我试图解决11个非线性联立方程与11个变量，其中所有变量都是积极的。
即使设置了所有解都应该为正的约束条件，我仍然得到否定的答案。例如，我设置f[6]= Pz，我理解为Pz〉=0，然而，最小化输出的Pz解是-1.38777878e-17。为什么会发生这种情况？我该怎么做才能解决这个问题（建议提出其他求解算法，如Nelder米德，为我的一套方程与这些约束是受欢迎的）。
我注意到，当我改变初始猜测时，解不再是负的，但我仍然认为，当我已经通过不等式约束指定所有解都是正的时，解不应该给出负的答案。

import numpy as np
from scipy.optimize import minimize
import sympy as sym
import math as math

###defining argument values
Cpz= 2 #argument
intensity=  0.04996599987353118 #argument
length = 10 #argument
Ax, Ay, Az, Bx, By, Bz, Cx, Cy, Cz, Dx, Dy, Dz = 0, length, 0, length, length,0,  0, 0, 0, length, 0, 0
ni, nm= 1.1, 1 #refrac index of optical grease and scint crystal #argument
thetacrit= sym.asin(nm/ni) #gives in radians
print (f"thetacrit in radians is {thetacrit}")

def my_fun(param):
    thetaC1= param[0]
    thetaC2= param[1]
    Cpx= param[2]
    Px=param[3]
    Pz=param[4]
    Cphorizdist= param[5]
    Cpvertdist= param[6]
    Cpdist3D= param[7]
    Chorizdist= param[8]
    Cvertdist= param[9]
    Cdist3D= param[10]
    f= np.zeros(11)
    
    f[0]= sym.asin((nm/ni)*sym.sin(thetaC2))- thetaC1
    f[1]= np.absolute(Px-Cpx)- Cphorizdist
    f[2]= np.absolute(Pz-Cpz)- Cpvertdist
    f[3]=( (Cphorizdist)**2 + (Cpvertdist)**2 )**(1/2)-Cpdist3D
    f[4]= np.absolute(Cpx-Cx)- Chorizdist
    f[5]= np.absolute(Cpz-Cz)-Cvertdist
    f[6]= (Chorizdist**2 + Cvertdist**2)**(1/2)- Cdist3D 
    f[7]= Cphorizdist/Cpdist3D-sym.sin(thetaC1)
    f[8]= Cpx/Cdist3D- sym.sin(thetaC2)
    f[9]= Cphorizdist/Cpvertdist- sym.tan(thetaC1)
    f[10]= 1/((Cpdist3D+Cdist3D)**2)-intensity
    
    return np.dot(f,f) #maybe add more

def my_cons(param):
    thetaC1= param[0]
    thetaC2= param[1]
    Cpx= param[2]
    Px=param[3]
    Pz=param[4]
    Cphorizdist= param[5]
    Cpvertdist= param[6]
    Cpdist3D= param[7]
    Chorizdist= param[8]
    Cvertdist= param[9]
    Cdist3D= param[10]

    f = np.zeros(13)
  
    #bc for solving method SLSQP, constraints dict type is ineq, these f[] become >=0
    #f[1] = -thetaC1

    
    f[0]= thetaC1
    f[1]= math.pi/2-thetaC1
    f[2]= thetaC2
    f[3]= math.pi/2-thetaC2
    f[4]= Cpx
    f[5]= Px
    f[6]= Pz
    f[7]= Cphorizdist
    f[8]= Cpvertdist
    f[9]= Cpdist3D
    f[10]= Chorizdist
    f[11] = Cvertdist
    f[12] = Cdist3D
    return f

cons = {'type' : 'ineq', 'fun': my_cons}
res = minimize(my_fun, (0.5,0.5,0.5,0.5,0.5,0.5,0.5,0.5,0.5,0.5,0.5), method='SLSQP',\
           constraints=cons,options= {"maxiter": 5000})
print(f"thetaC1,thetaC2,Cpx,Px,Pz,Cphorizdist,Cpvertdist,Cpdist3D,\
Chorizdist,Cvertdist,Cdist3D is {res}")