有解的非线性方程组上，西皮·解法失败了

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)

我尝试过改变fSolutions的参数，但没有效果。有人能帮我解决问题吗，这样求解者就能找到我正在寻找的解决方案？干杯。

Answer 1

我用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]

Answer 2

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

将函数更改为

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。