Scipy Fsolve在有解的非线性方程组上失败了

Question

我有一个由12个非线性方程组成的系统，如下所示：

import math
import numpy as np
# First: Defining some variables that have been calculated before (for the code to work)
n_s1=[8.71557427e-02, 9.96194698e-01, 6.12323400e-17]
n_s2=[5.23359562e-02, 9.98629535e-01, 6.12323400e-17]
P_sens1=[ 2.00000000e+00,  5.01223926e-01, -1.43937571e-17]
P_sens2=[ 2.00000000e+00,  4.20537892e-01, -2.17255041e-17]
# My 12 unknowns, initialized to a known solution: 
ep1_x,ep1_y,ep2_x,ep2_y,er1_x,er1_y,er2_x,er2_y,ep_x,ep_y,en_x,en_y = [0.0287630388891315, 
0.32876303888913166, 0.016591877224378986, 0.3165918772243792, 0.9961946980917455, 0.08715574274765804,
0.9986295347545739, 0.052335956242943855, 0.0287630388891315, 0.32876303888913166, -0.7071067811865475,
0.7071067811865476]
# My 12 equations: 
eq = 12*[0]
eq[0] = (ep_x*en_x+ep_y*en_y)/(en_x*n_s1[0]+en_y*n_s1[1])*n_s1[0]-ep1_x
eq[1] = (ep_x*en_x+ep_y*en_y)/(en_x*n_s1[0]+en_y*n_s1[1])*n_s1[1]-ep1_y
eq[2] = n_s1[0]-2*(en_x*n_s1[0]+en_y*n_s1[1])/np.linalg.norm(n_s1)*en_x-er1_x
eq[3] = n_s1[1]-2*(en_x*n_s1[0]+en_y*n_s1[1])/np.linalg.norm(n_s1)*en_y-er1_y
eq[4] = (P_sens1[0]-ep1_x)/math.sqrt(math.pow((P_sens1[0]-ep1_x),2)+math.pow((P_sens1[1]-ep1_y),2))-er1_x
eq[5] = (P_sens1[1]-ep1_y)/math.sqrt(math.pow((P_sens1[0]-ep1_x),2)+math.pow((P_sens1[1]-ep1_y),2))-er1_y
eq[6] = (ep_x*en_x+ep_y*en_y)/(en_x*n_s2[0]+en_y*n_s2[1])*n_s2[0]-ep2_x
eq[7] = (ep_x*en_x+ep_y*en_y)/(en_x*n_s2[0]+en_y*n_s2[1])*n_s2[1]-ep2_y
eq[8] = n_s2[0]-2*(en_x*n_s2[0]+en_y*n_s2[1])/np.linalg.norm(n_s2)*en_x-er2_x
eq[9] = n_s2[1]-2*(en_x*n_s2[0]+en_y*n_s2[1])/np.linalg.norm(n_s2)*en_y-er2_y
eq[10] = (P_sens2[0]-ep2_x)/math.sqrt(math.pow((P_sens2[0]-ep2_x),2)+math.pow((P_sens2[1]-ep2_y),2))-er2_x
eq[11] = (P_sens2[1]-ep2_y)/math.sqrt(math.pow((P_sens2[0]-ep2_x),2)+math.pow((P_sens2[1]-ep2_y),2))-er2_y

我知道系统有一个解决方案，因为在这种情况下，当我打印eq列表时，所有条目的计算结果都为零(即最小化(。

print(np.round(eq,9))
>> [-0. -0. -0. -0.  0. -0. -0. -0.  0.  0.  0. -0.]

然而，使用如下实现的scipyfsolve方法解决系统对我的不起作用

from scipy.optimize import fsolve
def equations(p):
ep1_x,ep1_y,ep2_x,ep2_y,er1_x,er1_y,er2_x,er2_y,ep_x,ep_y,en_x,en_y = p
return (eq)
print(fsolve(equations, 12*[0.5]))

我在读取时收到错误消息

>> RuntimeWarning: The iteration is not making good progress, as measured by the 
improvement from the last ten iterations.
warnings.warn(msg, RuntimeWarning)

我尝试过更改fsolve的参数，但没有效果。有人能帮我吗？这样求解器就能找到我想要的解决方案吗？干杯

Answer 1

您不计算点p的方程，因为您在equation函数之外定义了方程。

将功能更改为

def equations(p):
eq = np.zeros(p.shape[0])
ep1_x,ep1_y,ep2_x,ep2_y,er1_x,er1_y,er2_x,er2_y,ep_x,ep_y,en_x,en_y = p
eq[0] = (ep_x*en_x+ep_y*en_y)/(en_x*n_s1[0]+en_y*n_s1[1])*n_s1[0]-ep1_x
eq[1] = (ep_x*en_x+ep_y*en_y)/(en_x*n_s1[0]+en_y*n_s1[1])*n_s1[1]-ep1_y
eq[2] = n_s1[0]-2*(en_x*n_s1[0]+en_y*n_s1[1])/np.linalg.norm(n_s1)*en_x-er1_x
eq[3] = n_s1[1]-2*(en_x*n_s1[0]+en_y*n_s1[1])/np.linalg.norm(n_s1)*en_y-er1_y
eq[4] = (P_sens1[0]-ep1_x)/math.sqrt(math.pow((P_sens1[0]-ep1_x),2)+math.pow((P_sens1[1]-ep1_y),2))-er1_x
eq[5] = (P_sens1[1]-ep1_y)/math.sqrt(math.pow((P_sens1[0]-ep1_x),2)+math.pow((P_sens1[1]-ep1_y),2))-er1_y
eq[6] = (ep_x*en_x+ep_y*en_y)/(en_x*n_s2[0]+en_y*n_s2[1])*n_s2[0]-ep2_x
eq[7] = (ep_x*en_x+ep_y*en_y)/(en_x*n_s2[0]+en_y*n_s2[1])*n_s2[1]-ep2_y
eq[8] = n_s2[0]-2*(en_x*n_s2[0]+en_y*n_s2[1])/np.linalg.norm(n_s2)*en_x-er2_x
eq[9] = n_s2[1]-2*(en_x*n_s2[0]+en_y*n_s2[1])/np.linalg.norm(n_s2)*en_y-er2_y
eq[10] = (P_sens2[0]-ep2_x)/math.sqrt(math.pow((P_sens2[0]-ep2_x),2)+math.pow((P_sens2[1]-ep2_y),2))-er2_x
eq[11] = (P_sens2[1]-ep2_y)/math.sqrt(math.pow((P_sens2[0]-ep2_x),2)+math.pow((P_sens2[1]-ep2_y),2))-er2_y
return eq

并使用初始点(0.1, ..., 0.1)，即

sol = fsolve(equations, x0=0.1*np.ones(12))

给出了具有目标值的解决方案

array([-2.01876016e-12,  2.43202680e-11, -4.70620209e-11,  4.23225760e-11,
4.71587214e-11, -4.30741137e-11,  1.50569834e-12, -2.75213741e-11,
-2.21203056e-11, -5.79871359e-11,  2.21933583e-11,  5.72545206e-11])

请注意，警告消息仍然存在，只是表示在过去的十次迭代中没有取得好的进展。为了获得更好的收敛性，可以将雅可比算子传递给fsolve。

Answer 2

我使用python的scipy.optimize.leastsq解决了这个问题。详细的代码如下所示。此外，我可以将它用于我未来可能拥有的超定系统：

def f(x):
ep1_x,ep1_y,ep2_x,ep2_y,er1_x,er1_y,er2_x,er2_y,ep_x,ep_y,en_x,en_y = x
return np.asarray(((ep_x*en_x+ep_y*en_y)/(en_x*n_s1[0]+en_y*n_s1[1])*n_s1[0]-ep1_x,
(ep_x*en_x+ep_y*en_y)/(en_x*n_s1[0]+en_y*n_s1[1])*n_s1[1]-ep1_y,
n_s1[0]-2*(en_x*n_s1[0]+en_y*n_s1[1])/np.linalg.norm(n_s1)*en_x-er1_x,
n_s1[1]-2*(en_x*n_s1[0]+en_y*n_s1[1])/np.linalg.norm(n_s1)*en_y-er1_y,
(P_sens1[0]-ep1_x)/math.sqrt(math.pow((P_sens1[0]-ep1_x),2)+math.pow((P_sens1[1]-ep1_y),2))-er1_x,
(P_sens1[1]-ep1_y)/math.sqrt(math.pow((P_sens1[0]-ep1_x),2)+math.pow((P_sens1[1]-ep1_y),2))-er1_y,
(ep_x*en_x+ep_y*en_y)/(en_x*n_s2[0]+en_y*n_s2[1])*n_s2[0]-ep2_x,
(ep_x*en_x+ep_y*en_y)/(en_x*n_s2[0]+en_y*n_s2[1])*n_s2[1]-ep2_y,
n_s2[0]-2*(en_x*n_s2[0]+en_y*n_s2[1])/np.linalg.norm(n_s2)*en_x-er2_x,
n_s2[1]-2*(en_x*n_s2[0]+en_y*n_s2[1])/np.linalg.norm(n_s2)*en_y-er2_y,
(P_sens2[0]-ep2_x)/math.sqrt(math.pow((P_sens2[0]-ep2_x),2)+math.pow((P_sens2[1]-ep2_y),2))-er2_x,
(P_sens2[1]-ep2_y)/math.sqrt(math.pow((P_sens2[0]-ep2_x),2)+math.pow((P_sens2[1]-ep2_y),2))-er2_y))
def system(x,b):
return (f(x)-b)
b = np.zeros(12)
solution = optimize.leastsq(system, np.ones(12), args=b)[0]