如何使用具有非线性约束的 LMFIT 解决最小二乘问题?
How to solve a least square problem with LMFIT with non-linear constraints?
我手头有一个最小二乘问题:Ax=b
没有任何限制,可以像下面这样解决:
import numpy as np
from lmfit import Model, Parameters, Minimizer, report_fit
from numpy.linalg import lstsq
A = np.array([[ 143, -144, -343, 56, 98, 54]
,[-7632, 77, -277, 63, 999, 345 ]
,[ 55, -75, 74, 744, 765, 77]
,[-22 , 177, -28, 12, 54, 577]
,[-848 , -433 , 121, 54, 643, 88]
,[98, 75, 155, 87, 23, 54]
,[12, 56, 44, 44, 86, 46]
,[75, 22, 111, 965, 486, 345]])
b = [-114, -734, -270, 577, 676, 122, 245, 654]
b = np.reshape(b, (8, 1))
x, residuals, rnk, s = lstsq(A, b, rcond=-1)
print (" x =", x)
x = [[-0.04886372]
[-2.38958265]
[ 1.44719216]
[ 0.42944648]
[-1.23440929]
[ 1.98725466]]
我想在每对结果上添加一个重复约束,并使用 LMFIT 解决这个问题,因为定义约束似乎更方便。我试过了,但没用。约束是,例如,对于第一对结果,x1
必须为正数,并且 -x1 <= x2 <= x1
和下一对相同的模式。
def fcn2min(params, A, b):
x1 = params['x1']
x2 = params['x2']
x3 = params['x3']
x4 = params['x4']
x5 = params['x5']
x6 = params['x6']
x = np.array([x1, x2, x3, x4, x5, x6])
x = np.reshape(x, (6, 1))
return np.sum((np.dot(A, x)-b)**2)
params = Parameters()
params.add('x1', min=0)
params.add('xangle1', value=0.05, vary=True, min=(-np.pi/2)+(0.00000000001), max=(np.pi/2)-(0.00000000001))
params.add('x2', expr='(x1)*(sin(xangle1))')
params.add('x3', min=0)
params.add('xangle2', value=0.05, vary=True, min=(-np.pi/2)+(0.00000000001), max=(np.pi/2)-(0.00000000001))
params.add('x4', expr='(x3)*(sin(xangle2))')
params.add('x5', min=0)
params.add('xangle3', value=0.05, vary=True, min=(-np.pi/2)+(0.00000000001), max=(np.pi/2)-(0.00000000001))
params.add('x6', expr='(x5)*(sin(xangle3))')
mini = Minimizer(fcn2min, params, fcn_args=(A, b))
out = mini.leastsq()
print (" report_fit(result) =", report_fit(out.params))
首先,它不起作用,我很难定义和解决问题。其次,如果 x
有 100 或 1000 个元素,手动定义大规模问题的参数和约束将非常困难。
处理此问题的任何帮助表示赞赏。谢谢你。
说出你实际得到了什么总是有帮助的。说"it doesn't work"基本没用。您应该报告输出和收到的任何错误消息。包括对您认为应该得到什么的描述有时也会有所帮助。
几点:
- 您需要 NOT return
np.sum((np.dot(A, x)-b)**2)
和 return np.dot(A, x)-b
。拟合将进行平方和求和。
- 您应该为可变参数指定明确的初始值。
有了它,您将得到一个运行到完成并满足您的约束的拟合:
def fcn2min(params, A, b):
x1 = params['x1']
x2 = params['x2']
x3 = params['x3']
x4 = params['x4']
x5 = params['x5']
x6 = params['x6']
x = np.array([x1, x2, x3, x4, x5, x6]).reshape((6, 1))
return np.dot(A, x) - b
params = Parameters()
params.add('x1', 0.25, min=0)
params.add('xangle1', value=0.05, vary=True, min=(-np.pi/2)+(0.00000000001), max=(np.pi/2)-(0.00000000001))
params.add('x2', expr='(x1)*(sin(xangle1))')
params.add('x3', 0.25, min=0)
params.add('xangle2', value=0.05, vary=True, min=(-np.pi/2)+(0.00000000001), max=(np.pi/2)-(0.00000000001))
params.add('x4', expr='(x3)*(sin(xangle2))')
params.add('x5', 0.25, min=0)
params.add('xangle3', value=0.05, vary=True, min=(-np.pi/2)+(0.00000000001), max=(np.pi/2)-(0.00000000001))
params.add('x6', expr='(x5)*(sin(xangle3))')
mini = Minimizer(fcn2min, params, fcn_args=(A, b))
out = mini.leastsq()
report_fit(out, min_correl=0.5)
这将打印出:
[[Fit Statistics]]
# fitting method = leastsq
# function evals = 75
# data points = 8
# variables = 6
chi-square = 597869.757
reduced chi-square = 298934.879
Akaike info crit = 101.773493
Bayesian info crit = 102.250143
[[Variables]]
x1: 0.21393627 +/- 0.19714856 (92.15%) (init = 0.25)
xangle1: -0.01575861 +/- 9.59829799 (60908.27%) (init = 0.05)
x2: -0.00337120 +/- 2.05111043 (60842.16%) == '(x1)*(sin(xangle1))'
x3: 0.89232913 +/- 1.36874096 (153.39%) (init = 0.25)
xangle2: -0.86188487 +/- 2.50527077 (290.67%) (init = 0.05)
x4: -0.67734114 +/- 0.99201291 (146.46%) == '(x3)*(sin(xangle2))'
x5: 0.89482649 +/- 1.32577015 (148.16%) (init = 0.25)
xangle3: 1.57078737 +/- 398996.439 (25401047.11%) (init = 0.05)
x6: 0.89482649 +/- 4.38420893 (489.95%) == '(x5)*(sin(xangle3))'
[[Correlations]] (unreported correlations are < 0.500)
C(x1, x5) = 0.878
C(x5, xangle3) = 0.853
C(xangle1, xangle3) = 0.829
C(xangle1, x5) = 0.808
C(x1, xangle2) = -0.802
C(x3, xangle2) = 0.732
C(xangle2, x5) = -0.729
C(x1, xangle1) = 0.663
C(x1, xangle3) = 0.640
C(xangle2, xangle3) = -0.548
C(xangle1, xangle2) = -0.531
不确定性非常非常高,表明您的约束模型可能存在某些问题。
我手头有一个最小二乘问题:Ax=b
没有任何限制,可以像下面这样解决:
import numpy as np
from lmfit import Model, Parameters, Minimizer, report_fit
from numpy.linalg import lstsq
A = np.array([[ 143, -144, -343, 56, 98, 54]
,[-7632, 77, -277, 63, 999, 345 ]
,[ 55, -75, 74, 744, 765, 77]
,[-22 , 177, -28, 12, 54, 577]
,[-848 , -433 , 121, 54, 643, 88]
,[98, 75, 155, 87, 23, 54]
,[12, 56, 44, 44, 86, 46]
,[75, 22, 111, 965, 486, 345]])
b = [-114, -734, -270, 577, 676, 122, 245, 654]
b = np.reshape(b, (8, 1))
x, residuals, rnk, s = lstsq(A, b, rcond=-1)
print (" x =", x)
x = [[-0.04886372]
[-2.38958265]
[ 1.44719216]
[ 0.42944648]
[-1.23440929]
[ 1.98725466]]
我想在每对结果上添加一个重复约束,并使用 LMFIT 解决这个问题,因为定义约束似乎更方便。我试过了,但没用。约束是,例如,对于第一对结果,x1
必须为正数,并且 -x1 <= x2 <= x1
和下一对相同的模式。
def fcn2min(params, A, b):
x1 = params['x1']
x2 = params['x2']
x3 = params['x3']
x4 = params['x4']
x5 = params['x5']
x6 = params['x6']
x = np.array([x1, x2, x3, x4, x5, x6])
x = np.reshape(x, (6, 1))
return np.sum((np.dot(A, x)-b)**2)
params = Parameters()
params.add('x1', min=0)
params.add('xangle1', value=0.05, vary=True, min=(-np.pi/2)+(0.00000000001), max=(np.pi/2)-(0.00000000001))
params.add('x2', expr='(x1)*(sin(xangle1))')
params.add('x3', min=0)
params.add('xangle2', value=0.05, vary=True, min=(-np.pi/2)+(0.00000000001), max=(np.pi/2)-(0.00000000001))
params.add('x4', expr='(x3)*(sin(xangle2))')
params.add('x5', min=0)
params.add('xangle3', value=0.05, vary=True, min=(-np.pi/2)+(0.00000000001), max=(np.pi/2)-(0.00000000001))
params.add('x6', expr='(x5)*(sin(xangle3))')
mini = Minimizer(fcn2min, params, fcn_args=(A, b))
out = mini.leastsq()
print (" report_fit(result) =", report_fit(out.params))
首先,它不起作用,我很难定义和解决问题。其次,如果 x
有 100 或 1000 个元素,手动定义大规模问题的参数和约束将非常困难。
处理此问题的任何帮助表示赞赏。谢谢你。
说出你实际得到了什么总是有帮助的。说"it doesn't work"基本没用。您应该报告输出和收到的任何错误消息。包括对您认为应该得到什么的描述有时也会有所帮助。
几点:
- 您需要 NOT return
np.sum((np.dot(A, x)-b)**2)
和 returnnp.dot(A, x)-b
。拟合将进行平方和求和。 - 您应该为可变参数指定明确的初始值。
有了它,您将得到一个运行到完成并满足您的约束的拟合:
def fcn2min(params, A, b):
x1 = params['x1']
x2 = params['x2']
x3 = params['x3']
x4 = params['x4']
x5 = params['x5']
x6 = params['x6']
x = np.array([x1, x2, x3, x4, x5, x6]).reshape((6, 1))
return np.dot(A, x) - b
params = Parameters()
params.add('x1', 0.25, min=0)
params.add('xangle1', value=0.05, vary=True, min=(-np.pi/2)+(0.00000000001), max=(np.pi/2)-(0.00000000001))
params.add('x2', expr='(x1)*(sin(xangle1))')
params.add('x3', 0.25, min=0)
params.add('xangle2', value=0.05, vary=True, min=(-np.pi/2)+(0.00000000001), max=(np.pi/2)-(0.00000000001))
params.add('x4', expr='(x3)*(sin(xangle2))')
params.add('x5', 0.25, min=0)
params.add('xangle3', value=0.05, vary=True, min=(-np.pi/2)+(0.00000000001), max=(np.pi/2)-(0.00000000001))
params.add('x6', expr='(x5)*(sin(xangle3))')
mini = Minimizer(fcn2min, params, fcn_args=(A, b))
out = mini.leastsq()
report_fit(out, min_correl=0.5)
这将打印出:
[[Fit Statistics]]
# fitting method = leastsq
# function evals = 75
# data points = 8
# variables = 6
chi-square = 597869.757
reduced chi-square = 298934.879
Akaike info crit = 101.773493
Bayesian info crit = 102.250143
[[Variables]]
x1: 0.21393627 +/- 0.19714856 (92.15%) (init = 0.25)
xangle1: -0.01575861 +/- 9.59829799 (60908.27%) (init = 0.05)
x2: -0.00337120 +/- 2.05111043 (60842.16%) == '(x1)*(sin(xangle1))'
x3: 0.89232913 +/- 1.36874096 (153.39%) (init = 0.25)
xangle2: -0.86188487 +/- 2.50527077 (290.67%) (init = 0.05)
x4: -0.67734114 +/- 0.99201291 (146.46%) == '(x3)*(sin(xangle2))'
x5: 0.89482649 +/- 1.32577015 (148.16%) (init = 0.25)
xangle3: 1.57078737 +/- 398996.439 (25401047.11%) (init = 0.05)
x6: 0.89482649 +/- 4.38420893 (489.95%) == '(x5)*(sin(xangle3))'
[[Correlations]] (unreported correlations are < 0.500)
C(x1, x5) = 0.878
C(x5, xangle3) = 0.853
C(xangle1, xangle3) = 0.829
C(xangle1, x5) = 0.808
C(x1, xangle2) = -0.802
C(x3, xangle2) = 0.732
C(xangle2, x5) = -0.729
C(x1, xangle1) = 0.663
C(x1, xangle3) = 0.640
C(xangle2, xangle3) = -0.548
C(xangle1, xangle2) = -0.531
不确定性非常非常高,表明您的约束模型可能存在某些问题。