雅可比矩阵的伪逆问题

Question

我正在使用基于梯度的变异算法（Takahama，2006 年）实施 E-Constrained 差分进化。对于变异部分，我必须计算试验向量的每个约束函数的梯度，从而定义雅可比矩阵。我是通过有限差分来做的。我知道矩阵应该是 m x n，其中 m = 约束数，n = 维数（试验向量的长度）。试验向量的扰动将是 pseudo-inverse_of_Jacobian x delta_C。我不知道出了什么问题，但我收到消息 operands could not be broadcast with shapes (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

文章正文：

Answer 1

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

trial = trial-inv_Jc*delta_C

这将尝试逐个元素相乘（参见 this post）。我想你反而想做点积

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

或

trial = trial-inv_Jc.dot(delta_C)

Answer 2

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

雅可比矩阵的伪逆问题

Problems with pseudo-inverse of Jacobian matrix

python

differential-evolution