1

我正在使用基于梯度的突变算法（Takahama，2006）实现 E 约束差分进化。对于突变部分，我必须计算试验向量的每个约束函数的梯度，从而定义雅可比矩阵。我正在通过有限的差异来做到这一点。我理解矩阵应该是 mxn，其中 m = 约束数和 n = 维数（试验向量的长度）。试验向量上的扰动将是pseudo-inverse_of_Jacobian x delta_C。我不知道出了什么问题，但我收到消息操作数无法与形状 (4,5) (5,1) 一起广播。你能帮忙吗？这个乘法不应该产生一个 4x1 数组吗？谢谢。

trial = np.zeros(4)
inv_Jc = np.linalg.pinv(Jacobian(trial))
delta_C = np.array(calculate_C(trial))
delta_C = delta_C.reshape(len(delta_C),1)
trial = trial-inv_Jc*delta_C

您可以在下面找到我要解决的示例问题的相关定义。

import numpy as np

def fobj(x): 
    return 3*x[0]+0.000001*x[0]**3+2*x[1]+(0.000002/3)*x[1]**3

def constraints(x): 
    g = [-x[3]+x[2]-0.55,
         -x[2]+x[3]-0.55]
    h = [1000*np.sin(-x[2]-0.25)+1000*np.sin(-x[3]-0.25)+894.8-x[0],
         1000*np.sin(x[2]-0.25)+1000*np.sin(x[2]-x[3]-0.25)+894.8-x[1],
         1000*np.sin(x[3]-0.25)+1000*np.sin(x[3]-x[2]-0.25)+1294.8]
    
    return np.array([max(0,g[i]) for i in range(len(g))]), np.array(h)
     
def calculate_C(v):
    g, h = constraints(v)
    return np.array([*g,*h])

def Jacobian(v,delta=0.01):
    C = calculate_C(v)
    Jc = np.zeros(shape=(len(C),len(v)))

    for i in range(len(delta_C)):
        for j in range(len(v)):
            aux = np.copy(v)
            aux[j] = v[j]+delta
            C_plus = calculate_C(aux)[i]
            aux[j] = v[j]-delta
            C_minus = calculate_C(aux)[i]
            Jc[i,j] = ((C_plus-C_minus)/(2*delta))

    return Jc

文章中的文字：

基于梯度的突变

4

2 回答 2

1

我认为您的问题是您*在这一行中使用了运算符

trial = trial-inv_Jc*delta_C

这将尝试逐个元素相乘（参见这篇文章）。我认为您反而想采用点积

trial = trial-np.dot(inv_Jc, delta_C)

或者

trial = trial-inv_Jc.dot(delta_C)
于 2021-06-22T23:04:07.413 回答
0

@falafelocelot 向我暗示了正确的方向。验证数组的形状时，我注意到我正在使用定义为 的试验向量测试函数trial = np.zeros(4)，其形状为 (4,)。正如预期的那样，使用np.dot()提供了一个 (4,1) 数组，但该操作trial-np.dot(a,b)提供了一个 (4,4) 数组。当我将试验向量重新整形为 (4,1) 时，一切都按预期工作。

于 2021-06-23T12:10:39.293 回答