1

感谢您的关注，我已经听取了建议并尽可能简化了代码以隔离故障。

最小化中使用的函数找不到我的频率或幅度断点。它似乎与将 getpaz() def 及其最小化函数嵌入另一个定义中有关。当我在主循环中使用最小化代码运行 Getpaz() 函数时，它会运行。当我创建定义（例如通道处理 def）时，它失败了。

应该运行以下代码，在 Jupyter Notebooks 以及在命令行中编译的 python 3.x 中重复该问题。问题是，最小化函数找不到我引用的变量“响应”或“频率”，即使它们明显存在。我不明白为什么最小化功能会失败。这个精彩的节目让我倒退了三个星期。我无计可施，试图弄清楚。

如果我采用定义 paztest_fixed() 并将其代码剥离并将其放入可执行循环中，它就会运行。有没有人了解导致它失败的原因以及我如何解决它？我真的需要一个定义，以便我可以为整个历史地震台站目录处理多个通道！

import numpy as np
import scipy.signal
import scipy.optimize


# Daniel Burk, Michigan State University


def pazto_freq_resp(freqs, zeros, poles, scale_fac):
    b, a = scipy.signal.ltisys.zpk2tf(zeros, poles, scale_fac)
    if not isinstance(a, np.ndarray) and a == 1.0:
        a = [1.0]
    return scipy.signal.freqs(b, a, freqs * 2 * np.pi)[1] 
            # list element 0 is frequencies
            # list element 1 is the complex amplitudes

def phasecalc(testresponse):  # Bring in a list of complex numbers and return the angle between 90 and 270 degrees
    testphase = []
    for t in testresponse:
        tp = np.arctan2(t.imag , t.real) * 180. / np.pi
        if tp > 90.:
            tp = tp-360.
        testphase.append(tp - 90.0) # adjust phase to better match what is seen in the real world calibrations
    return(testphase)


# This is the function definition that is used in the scipy.minimize function.
def minimize(_var):     # Uses data found in frequencies, and in response. 
                        # Make sure phase and response tables use the same subset of frequencies.
    p1r, p1i, p3r, p4r, p5r,z1r,z2r,z3r, scale_fac = _var
    new_resp = pazto_freq_resp(
        freqs=frequencies,
        zeros=np.array([0.0 + 0.0 * 1j,
                        0.0 + 0.0 * 1j,
                        z1r + 0.0 * 1j,
                        z2r + 0.0 * 1j,
                        z3r + 0.0 * 1j], dtype=np.complex128),                        
        poles=np.array([p1r + p1i * 1j,
                        p1r - p1i * 1j,
                        p3r + 0.0 * 1j,
                        p4r + 0.0 * 1j,
                        p5r + 0.0 * 1j], dtype=np.complex128),
        scale_fac=scale_fac)
    return ((np.abs(new_resp) - np.abs(response)) ** 2).sum()



def getpaz(frequencies,response,Phasedeg):

    evaluation = 1.0E+09 # For evaluating how close the solution is to the original curve
    paz = [] # The new poles and zeros for the channels

    for z in range(0,32): # iterate 32 times to find the solution that best describes the phase response.
        initial_x=[]
        X0=np.random.random(9)
        #                                Using the minimize function, find the poles & zeros solution that best describes
        #                                the instrument response as found in responses, on frequencies breakpoint "frequencies"
        out = scipy.optimize.minimize(
            fun=minimize,
            method="BFGS",
            x0=X0,
            options={"eps": 1e-10, "maxiter": 1e8})

        x = out.x
        new_poles = np.array([-abs(x[0]) + abs(x[1]) * 1j,
                              -abs(x[0]) - abs(x[1]) * 1j, 
                              -abs(x[2]) + 0.0 * 1j,
                              -abs(x[3]) + 0.0 * 1j,
                              -abs(x[4]) + 0.0 * 1j], 
                              dtype=np.complex128)    

        new_zeros = np.array([ 0.0 + 0.0 * 1j,
                               0.0 + 0.0 * 1j,
                              x[5] + 0.0 * 1j,
                              x[6] + 0.0 * 1j,
                              x[7] + 0.0 * 1j], dtype=np.complex128)
        new_scale_fac = x[8]
        #              Create the response curve that results from this theoretical new poles and zeroes solution
        inverted_response = pazto_freq_resp(freqs=frequencies, zeros=new_zeros, poles=new_poles,scale_fac=new_scale_fac)    
        inphase = phasecalc(inverted_response) # phase from inverted response, listed in degrees
        curvefit = np.sqrt(((np.array(Phasedeg) - np.array(inphase))**2).mean()) # rmse
        if (curvefit) < evaluation:
            final_iteration = z
            best_poles=new_poles
            best_zeros=new_zeros
            best_scale_fac=new_scale_fac
            print(f'Iteration # {z}: Phase misfit drops to {curvefit}')
            evaluation = curvefit
    return(best_poles,best_zeros,best_scale_fac,evaluation,z)

def paztest_fixed():
    #################################### test with def #################################
    Component = 'LM.NE8K.MHZ'
    Caldate   = '05/15/2019'
    #              Frequency breakpoints within the passband of the seismometer to simulate
    frequencies = np.array([0.05, 0.0571, 0.0667, 0.080, 0.111, 0.133, 0.167,    \
                   0.222, 0.250, 0.286, 0.333, 0.400, 0.500, 0.526,     \
                   0.555, 0.588, 0.625, 0.666, 0.714, 0.770, 0.833,     \
                   0.910, 1.000, 1.111, 1.250, 1.429, 1.667, 2.000,     \
                   3.000, 4.000, 5.000, 8.000])

    #               here are the gain values for the seismometer at the above frequencies
    response = np.array([3.00, 4.48, 7.11, 12.28, 32.81, 56.56, 110.00, 258.43,      \
            366.09, 542.60, 852.84, 1451.12, 2764.14, 3201.37, 3734.65, \
            4390.91, 5205.66, 6225.33, 7508.61, 9123.91, 11134.60,      \
            13556.01, 16267.16, 18911.45, 20951.61, 22004.53, 22120.93, \
            21630.53, 20132.44, 18990.64, 17947.77, 15053.22])

    #           phase delay of the light beam vs. ground motion in degrees at above frequencies
    Phasedeg = [-6.660, -7.714, -8.880, -10.800, -14.800, -18.000, -22.200, \
                -30.000, -33.300, -37.543, -44.400, -52.560, -66.600,       \
                -69.158, -73.800, -78.353, -83.475, -89.520, -96.429,       \
                -104.677, -114.600, -126.654, -140.760, -157.600, -175.500, \
                -193.886, -211.560, -228.024, -254.865, -269.640, -280.080, \
                -300.528]

    best_poles,best_zeros,best_scale_fac,evaluation,final_iteration = getpaz(frequencies,response,Phasedeg)

    print("\n========================================")
    print(f"Inverted values for {Component} on caldate of {Caldate}:")
    print("Poles =\n", best_poles)
    print("Zeros =\n", best_zeros)
    print("Scale_fac = ", best_scale_fac)
    print(f"Evaluated misfit of phase = {evaluation} on iteration # {final_iteration}\n")
    print("========================================")


##############################  RUN CODE AS A DEF ######################

paztest_fixed()
4

1 回答 1

0

从 scipy 文档（https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.minimize.html），我们看到

scipy.optimize.minimize(fun, x0, args=(), method=None, jac=None, hess=None, hessp=None, bounds=None, constraints=(), tol=None, callback=None, options=无）[来源]

和

args ：元组，可选。传递给目标函数及其导数（fun、jac 和 hess 函数）的额外参数。

您的最小化功能需要以某种方式知道什么frequencies和response是什么。修改您的最小化函数，使其接受这两个参数，然后在您的调用中添加一个参数minimize以包含这些参数。它应该看起来像这样：

  1. 改变minimize定义：
def minimize(_var, frequencies, response):
    # Rest of your code
  1. 将这些参数添加到您的调用中以最小化（在getpaz函数内部）
out = scipy.optimize.minimize(
            fun=minimize,
            args = (frequencies, response), # <--- This is what I added
            method="BFGS",
            x0=X0,
            options={"eps": 1e-10, "maxiter": 1e8})

这些变化给了我这个输出。

Iteration # 0: Phase misfit drops to 1.224353466864608

========================================
Inverted values for LM.NE8K.MHZ on caldate of 05/15/2019:
Poles =
 [ -3.13071301+6.04747984j  -3.13071301-6.04747984j
  -4.50400442+0.j         -32.91191804+0.j
 -52.80201143+0.j        ]
Zeros =
 [ 0.00000000e+00+0.j  0.00000000e+00+0.j -3.24109312e+01+0.j
  1.19677764e-04+0.j  6.81647031e+02+0.j]
Scale_fac =  1602.3460317044198
Evaluated misfit of phase = 1.224353466864608 on iteration # 31

========================================
于 2019-05-15T16:55:54.330 回答