I am trying to solve 11 non-linear simultaneous equations with 11 variables, where all variables are positive.
Even after setting the constraints that all solutions should be positive, I am still getting negative answers. For example I set f[6]= Pz which I understand to mean Pz>=0, however the solution for Pz from the minimize output is -1.38777878e-17. Why is this happening? What can I do to resolve this (suggestions proposing other solving algorithms eg. Nelder Mead that work for my set of eqns with these constraints are welcome).
I noticed the solutions are no longer negative when I change the initial guess, but I still think the solutions should not be giving negative answers when I have already specified all solutions to be positive through inequality constraints.
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}")
Your constraints are actually just simply bounds on the optimization variables, so I'd suggest to pass them as bounds instead of using the more general constraints:
x0 = 0.5*np.ones(12)
bounds = [(0, None) for _ in range(x0.size)]
bounds[0] = (0, np.pi/2) # thetaC1
bounds[1] = (0, np.pi/2) # thetaC2
cons = {'type' : 'ineq', 'fun': my_cons}
res = minimize(my_fun, x0=x0, bounds=bounds, method='SLSQP',\
constraints=cons,options= {"maxiter": 5000})
That being said, it's also worth mentioning that
you should try do avoid absolute values whenever possible. Your objective function is not continuously differentiable due to the absolute values. Consequently, you're violating the mathematical assumptions of the SLSQP algorithm, which can lead to really odd behaviour once the algorithm reaches a point of discontinuity (i.e. a point x where the absolute value is zero).
there's no real need to mix the sympy, numpy and math modules. Numpy already provides all the functions and constants you're using (np.arcins, np.sin, np.pi etc).