Scipy.Optimize中IndexError排查：数组索引过多错误

Hey, let's break down why you're hitting that IndexError and fix it step by step:

What's Causing the Error? #

There are two core issues here:

1. Your Constraint Function Isn't Receiving Required Arguments #

When you defined your equality constraint, you forgot to pass the args parameter to create_const:

con1 = {'type': 'eq', 'fun': create_const}

scipy.optimize.minimize doesn't automatically share the args from your target function with constraints. So when the optimizer calls create_const, it only sends the optimization variable dy—your *args ends up empty. That makes arg = np.asarray(args) an empty array, and trying to access arg[0, i] throws the index error.

2. Argument Parsing Logic is Broken #

When you tested create_const(dy0, arg), you passed the entire arg tuple as a single argument. This meant *args became (arg,), so np.asarray(args) was a 2D array (shape (1, N)), and arg[0, i] worked. But when you correctly pass args to the constraint, *args will be all the individual elements of arg (like type_x[0], type_x[1], ..., delta2[2]), making np.asarray(args) a 1D array. Trying to use 2D indexing ([0, i]) on a 1D array will fail.

How to Fix It #

Let's resolve both issues:

Step 1: Pass Arguments to the Constraint #

Update your constraint dictionary to include the args parameter, wrapping arg in a tuple so it's passed as a single argument to create_const:

con1 = {'type': 'eq', 'fun': create_const, 'args': (arg,)}

Step 2: Fix Argument Parsing in create_const #

Adjust the function to accept the entire arg tuple as a single parameter, then index it correctly (since it's now a 1D array):

def create_const(dy, arg):  # Remove *args, take arg as a single parameter
    arg = np.asarray(arg)
    n = np.shape(dy)[0]
    dx = np.zeros((n, 2))
    type_x = np.zeros(n, dtype=object)  # Use object dtype to hold string values
    delta1 = np.zeros(n)
    delta2 = np.zeros(n)
    gamma = np.zeros((n, n))
    
    for i in range(n):
        # arg[i] corresponds to type_x[i]
        a, b = bndr(arg[i])
        # delta1 starts at index n+1 (after n type_x elements + 1 dP value)
        delta1[i] = arg[n + 1 + i]
        # delta2 starts at index 2n+1
        delta2[i] = arg[2*n + 1 + i]
        dx[i, 0] = (b - a) * dy[i]
    
    # Note: dx[:,1] remains all zeros here—double-check if this is intentional for your logic
    gamma = GammaApprox(delta1, delta2, dx[:, 1], dx[:, 0])
    d = np.dot(delta2, dx[:, 0])
    g = np.dot(dx[:, 0], gamma)
    g = np.dot(g, dx[:, 0])
    dP = float(arg[n])  # dP is at index n (right after the n type_x elements)
    
    return d + 0.5 * g - dP

Step 3: Keep Your Test Call Valid #

Since we adjusted create_const to take arg as a single parameter, your existing test line works as-is:

testconst = create_const(dy0, arg)

Quick Side Check #

You might want to verify the dx[:,1] assignment in create_const—right now it's initialized to zero but never updated. That could lead to unexpected results in GammaApprox, so make sure that's intentional for your use case.

Full Modified Code #

Here's the complete code with all fixes applied:

from scipy.optimize import minimize
import numpy as np

def totaldist(dy):
    n = np.shape(dy)[0]
    temp = 0
    for i in range(n):
        temp += dy[i] ** 2
    return -0.5 * temp

def create_bond(dy_max):
    n = np.shape(dy_max)[0]
    bond = np.zeros((n, 2))
    for i in range(n):
        bond[i, :] = [0, dy_max[i]]
    tot = tuple([tuple(row) for row in bond])
    return tot

def create_const(dy, arg):
    arg = np.asarray(arg)
    n = np.shape(dy)[0]
    dx = np.zeros((n, 2))
    type_x = np.zeros(n, dtype=object)
    delta1 = np.zeros(n)
    delta2 = np.zeros(n)
    gamma = np.zeros((n, n))
    
    for i in range(n):
        a, b = bndr(arg[i])
        delta1[i] = arg[n + 1 + i]
        delta2[i] = arg[2*n + 1 + i]
        dx[i, 0] = (b - a) * dy[i]
    
    gamma = GammaApprox(delta1, delta2, dx[:, 1], dx[:, 0])
    d = np.dot(delta2, dx[:, 0])
    g = np.dot(dx[:, 0], gamma)
    g = np.dot(g, dx[:, 0])
    dP = float(arg[n])
    
    return d + 0.5 * g - dP

def GammaApprox(delta1, delta2, x1, x2):
    n = np.shape(delta1)[0]
    gamma = np.zeros((n, n))
    for i in range(n):
        if x2[i] == x1[i]:
            gamma[i, i] = 0
        else:
            gamma[i, i] = (delta2[i] - delta1[i]) / (x2[i] - x1[i])
    return gamma

def GetNewPoint(x1, x2, delta1, delta2, type_x, P):
    n = np.shape(delta1)[0]
    dmax = np.zeros(n)
    dy0 = np.zeros(n)
    
    for i in range(n):
        a, b = bndr(type_x[i])
        if x2[i] > x1[i]:
            dmax[i] = (x2[i] - x1[i])/(b - a)
            dy0[i] = 1 / (b - a) * (x2[i] - x1[i]) / 2
        else:
            dmax[i] = (x1[i] - x2[i])/(b - a)
            dy0[i] = 1 / (b - a) * (x1[i] - x2[i]) / 2
    
    bond = create_bond(dmax)
    arg = ()
    for i in range(n):
        arg = arg + (type_x[i],)
    arg = arg + (abs(P[0] - P[1]), )
    for i in range(n):
        arg = arg + (delta1[i], )
    for i in range(n):
        arg = arg + (delta2[i], )
    
    testconst = create_const(dy0, arg)
    con1 = {'type': 'eq', 'fun': create_const, 'args': (arg,)}
    cons = ([con1, ])
    
    solution = minimize(totaldist, dy0, args=arg, method='SLSQP', bounds=bond, constraints=cons, options={'disp': True})
    x = solution.x
    print(x)
    return x

def bndr(type_x):
    if type_x == 'normal':
        x_0 = -5
        x_f = 1.5
    if type_x == 'lognorm':
        x_0 = 0.0001
        x_f = 5
    if type_x == 'chisquare':
        x_0 = 0.0001
        x_f = (0.8 * (10 ** .5))
    return x_0, x_f

def test():
    x1 = np.array([0.0001, 0.0001, -5])
    x2 = np.array([1.6673, 0.84334, -5])
    delta1 = np.array([0, 0, 0])
    delta2 = np.array([2.44E-7, 2.41E-6, 4.07E-7])
    type_x = np.array(['lognorm', 'chisquare', 'normal'])
    P = (0, 6.54E-8)
    f = GetNewPoint(x1, x2, delta1, delta2, type_x, P)
    return f

test()

内容的提问来源于stack exchange，提问作者Benji_90210