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

如何解决Scipy Fsolve 在具有解的非线性方程组上失败

我有一个由 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.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))-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))-er2_y

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

print(np.round(eq,9))

>> [-0. -0. -0. -0.  0. -0. -0. -0.  0.  0.  0. -0.]

但是，使用如下实现的 scipy fsolve 方法解决系统对我不起作用

from scipy.optimize import fsolve

def equations(p):
    ep1_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 的参数但没有任何效果。任何人都可以帮助我，以便求解器可以找到我正在寻找的解决方案吗？干杯。

解决方法

您不计算点 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))-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))-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。

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

def f(x):
    ep1_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))-er1_x,(P_sens1[1]-ep1_y)/math.sqrt(math.pow((P_sens1[0]-ep1_x),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))-er2_x,(P_sens2[1]-ep2_y)/math.sqrt(math.pow((P_sens2[0]-ep2_x),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]