1

我正在尝试优化我编写的代码来计算方程组的最小二乘并返回未知数的最佳值：（a1,a2,a3,z1,z2并且pottemp是zlevels已知的）。

方程组如下（https://imgur.com/LQqeOgA）：

方程组

我编写了以下代码，它就像一个魅力，但它不是很有效，你可以看到由于 for 循环的数量，如果我增加 and 的数量hstep，astep我的数组变得非常大并且需要很多时间从 计算最小二乘leastsquared。

def leastsquared(zlevels,pottemp,z1,z2,a1,a2,a3):
    dtot=0.0
    for zi in range(len(zlevels)): #zvalue of obs points 
        if zlevels[zi]<=z1:
            F=np.tan(a1)*zlevels[zi]+pottemp[0]
        elif zlevels[zi]<=z2:
            F=(np.tan(a1)*z1+pottemp[0]) + np.tan(a2)*(zlevels[zi]-z1)
        else: #zlevels[zi]<=H
            F=((np.tan(a1)*z1+pottemp[0]) + np.tan(a2)*(z2-z1) ) + np.tan(a3)*(zlevels[zi]-z2)
        d=(pottemp[zi]-F)**2.0
        dtot+=d
    dtot=dtot/len(zlevels)
    #print("dtot is" + str(dtot))
    return dtot

def optm(hstep,astep,zlevels,pottemp):

    sqdist=np.inf
    for z1 in np.linspace(0,zlevels[-1],hstep):
        for z2 in np.linspace(0,zlevels[-1],hstep):
            for a1 in np.linspace(0,0.1,astep): # max angle
                for a2 in np.linspace(0,0.01,astep):
                    for a3 in np.linspace(0,0.01,astep):
                        sqdist_new=leastsquared(zlevels,pottemp,z1,z2,a1,a2,a3)
                        if sqdist_new<sqdist:
                            sqdist=sqdist_new
                            optimalsetting=[z1,z2,a1,a2,a3]

    return optimalsetting

代码本身非常简单但效率低下。我试图使用 实现此代码scipy.optimize.minimize，但我无法运行它。有谁知道如何优化它？也许，scipy.optimize.minimize这将是最简单的方法，但我不确定如何使我的代码适应它。

4

2 回答 2

1

肯定有一种更有效的方法来处理这个问题，因为对分段线性问题进行了很好的研究，但以下是scipy.optimize根据要求完成的。scipy 也无法处理动态约束，例如要求z1<z2，所以如果没有强制执行，我交换了求解器中变量的顺序。通过删除这些代码行并让数据自己说话，优化可能会更加稳定。

import numpy as np
import scipy.optimize

def f(zlevels, t0, a1, a2, a3, z1, z2, H):
    """
    Function to evaluate f, given a set of variables
    """
    # Check that no samples are out of bounds
    if any(zlevels < 0) or any(zlevels > H):
        quit("Values of zlevels out of bounds")

    pottemp = np.zeros(zlevels.shape)
    mask1 = zlevels <= z1
    mask2 = (z1 < zlevels) & (zlevels <= z2)
    mask3 = z2 < zlevels

    z1_arr1 = np.ones(sum(mask2))*z1
    z1_arr2 = np.ones(sum(mask3))*z1
    z2_arr = np.ones(sum(mask3))*z2

    pottemp[mask1] = zlevels[mask1] * a1 + t0
    pottemp[mask2] = z1_arr1 * a1 + t0 + a2 * (zlevels[mask2] - z1)
    pottemp[mask3] = z1_arr2 * a1 + t0 + a2 * (z2_arr - z1) + a3 * (zlevels[mask3] - z2)

    return pottemp

def obj_fun(x, zlevels, pottemp, t0, H):
    """
    Least-squares objective function for the problem
    """
    a1, a2, a3, z1, z2 = x
    # Swap order if z1 is larger than z2
    if z1 > z2:
        z1, z2 = z2, z1

    pottemp_pred = f(zlevels, t0, a1, a2, a3, z1, z2, H)
    return sum((pottemp_pred-pottemp)**2) 

def optimize(zlevels, pottemp, t0, H):
    """
    Optimize a1, a2, a3, z1, z2
    """

    starting_guess = [0,0,0,0.5,2.5]
    res = scipy.optimize.minimize(obj_fun, starting_guess, args=(zlevels, pottemp, t0, H), \
            bounds=[(None,None),(None,None),(None,None),(0,H),(0,H)])
    x = res.x
    # Swap order if z1 is larger than z2
    if x[-2] > x[-1]:
        x[-2], x[-1] = x[-1], x[-2]
    return x

# Make the true model for testing
t0 = 0.5
a1 = -1
a2 = +1
a3 = -2
z1 = 1
z2 = 2
H = 3

# Generate some variables
n = 1000
zlevels = np.random.random(n)*3
# Get values and add some random noise
pottemp = f(zlevels, t0, a1, a2, a3, z1, z2, H) + np.random.normal(0,0.1, size=n)

print("Correct values:", a1, a2, a3, z1, z2)
a1, a2, a3, z1, z2 = optimize(zlevels, pottemp, t0, H)
print("Fitted values:", a1, a2, a3, z1, z2)
于 2019-09-03T11:02:07.860 回答
0

我不太确定实际变量及其值。但下面是我尝试过的东西：

import numpy as np
from scipy.optimize import minimize

def leastsquared(zlevels,pottemp,z1,z2,a1,a2,a3):
    dtot=0.0
    for zi in range(len(zlevels)): #zvalue of obs points 
        if zlevels[zi]<=z1:
            F = np.tan(a1)*zlevels[zi]+pottemp[0]
        elif zlevels[zi]<=z2:
            F = (np.tan(a1)*z1+pottemp[0]) + np.tan(a2)*(zlevels[zi]-z1)
        else: #zlevels[zi]<=H
            F=((np.tan(a1)*z1+pottemp[0]) + np.tan(a2)*(z2-z1) ) + np.tan(a3)*(zlevels[zi]-z2)
        d=(pottemp[zi]-F)**2.0
        dtot+=d
    dtot=dtot/len(zlevels)
    #print("dtot is" + str(dtot))
    return dtot

def optm(hstep,astep,zlevels,pottemp):
    steps = np.linspace(0,0.01,astep)
    x0 = np.array([steps, steps, steps]) #a1, a2, a3

    def objective(x):
        return leastsquared(zlevels,pottemp, z1, z2, x[0], x[1], x[2])

    sqdist=np.inf
    for z1 in np.linspace(0,zlevels[-1],hstep):
        for z2 in np.linspace(0,zlevels[-1],hstep):
            sol = minimize(objective, x0, method='SLSQP', options={'disp':True})
            optimalsetting=[z1,z2,sol.x[0],sol.x[1],sol.x[2]]
    return optimalsetting


# test
# some random initialization
astep = 2
steps = 3


optm(5, 2, [1, 2], [1, 2])

您可能想在此处更新目标函数和变量类型。希望这能让你开始。

于 2019-09-02T11:53:24.993 回答