python-3.x - 尝试使用 GEKKO OPTIMIZER 时出现 "No loop matching the specified signature and casting was found for ufunc solve"

Question

我正在尝试使用 GEKKO 优化我的 FEA 桁架解决方案。所以我将所有 FEA 解决方案定义为一个函数，并尝试将力作为变量。所以目标是通过一个约束(所有力的总和 =-100000)最小化力。但是在我运行我的代码后我得到一个错误我无法修复它谷歌搜索也没有帮助(我是 python 的新手)。如果有人能告诉我问题出在哪里，我将不胜感激。提前致谢。

尝试: 来自 gekko 进口 GEKKO 除了: # pip 安装 gekko 导入点 pip.main(['安装','gekko']) 从 gekko 导入 GEKKO

# from __future__ import division
import numpy as np
import matplotlib.pyplot as plt

    # numElem=11
    # numNodes=6


def calculate_truss_force(force):
    # defined coordinate system
    x_axis = np.array([1, 0])
    y_axis = np.array([0, 1])

    # elements coordinates
    elemNodes = np.array([[0, 1], [0, 2], [1, 2], [1, 3],
                          [0, 3], [2, 3], [2, 5], [3, 4], [3, 5], [2, 4], [4, 5]])
    # nodes coordinates
    nodeCords = np.array([
        [0.0, 0.0], [0.0, 100.0],
        [100.0, 0.0], [100.0, 100.0],
        [200.0, 0.0], [200.0, 100.0]])


    modE = 200000
    Area = 200
    # assembling the model

    numElem = elemNodes.shape[0]
    numNodes = nodeCords.shape[0]



    xx = nodeCords[:, 0]
    yy = nodeCords[:, 1]

    EA = modE * Area
    tdof = 2 * numNodes  # total number of degrees of freedom
    disps = np.zeros((tdof, 1))
#    force = np.zeros((tdof, 1))
    sigma = np.zeros((numElem, 1))
    stiffness = np.zeros((tdof, tdof))
    np.set_printoptions(precision=3)

    # applying the load

#    force[3] = -20000.0
#    force[7] = -50000.0
#    force[11] = -30000.0


    # defined boundary
    presDof = np.array([0, 1, 9])
    for e in range(numElem):
        indice = elemNodes[e, :]
        elemDof = np.array([indice[0] * 2, indice[0] * 2 + 1, indice[1] * 2, indice[1] * 2 + 1]) #dof corresponding to the global stifnessmatrix
        xa = xx[indice[1]] - xx[indice[0]]   #length of the elemnts in x dir
        ya = yy[indice[1]] - yy[indice[0]]   #lenth of the elemnt in y dir
        len_elem = np.sqrt(xa * xa + ya * ya)
        c = xa / len_elem
        s = ya / len_elem
        k1 = (EA / len_elem) * np.array([[c * c, c * s, -c * c, -c * s],
                                         [c * s, s * s, -c * s, -s * s],
                                         [-c * c, -c * s, c * c, c * s],
                                         [-c * s, -s * s, c * s, s * s]])
        stiffness[np.ix_(elemDof, elemDof)] += k1    # add elem K to the global K

    actDof = np.setdiff1d(np.arange(tdof), presDof)

    disp1 = np.linalg.solve(stiffness[np.ix_(actDof, actDof)], force[np.ix_(actDof)])

    disps[np.ix_(actDof)] = disp1

    # stresses at elements

    for e in range(numElem):
        indice = elemNodes[e, :]
        elemDof = np.array([indice[0] * 2, indice[0] * 2 + 1, indice[1] * 2, indice[1] * 2 + 1])
        xa = xx[indice[1]] - xx[indice[0]]
        ya = yy[indice[1]] - yy[indice[0]]
        len_elem = np.sqrt(xa * xa + ya * ya)
        c = xa / len_elem
        s = ya / len_elem
        sigma[e] = (modE / len_elem) * np.dot(np.array([-c, -s, c, s]), disps[np.ix_(elemDof)])

    #print (disps)
    #print (sigma)
    return np.max(np.abs(sigma))
m = GEKKO()

# Define variables
A = m.Array(m.Var, (12))

# initial guess
ig = [0, 0, 0, -20000, 0, 0, 0, -50000, 0, 0, 0, -30000]

#  bounds
for i, Ai in enumerate(A):
    Ai.value = ig[i]
    Ai.lower = ig[i] * 0.95
    Ai.upper = ig[i] * 1.05
m.Equation(np.sum(A) == -100000)
m.Obj(calculate_truss_force(A.reshape(12, 1)))
m.solve()
print(A.reshape(12, 1))
print(calculate_truss_force(np.array(ig).reshape(12, 1)))

我切换到SCipy，代码如下

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import minimize

    # numElem=11
    # numNodes=6

def calculate_truss_force(Area_elem):
#    print("calc_truss")
    # defined coordinate system
    x_axis = np.array([1, 0])
    y_axis = np.array([0, 1])

    # elements coordinates
    elemNodes = np.array([[0, 1], [0, 2], [1, 2], [1, 3],
                          [0, 3], [2, 3], [2, 5], [3, 4], [3, 5], [2, 4], [4, 5]])
    # nodes coordinates
    nodeCords = np.array([
        [0.0, 0.0], [0.0, 100.0],
        [100.0, 0.0], [100.0, 100.0],
        [200.0, 0.0], [200.0, 100.0]])

    modE = 200000
    Area = 200
    # assembling the model

    numElem = elemNodes.shape[0]
    numNodes = nodeCords.shape[0]

    xx = nodeCords[:, 0]
    yy = nodeCords[:, 1]

    EA = modE * Area
    tdof = 2 * numNodes  # total number of degrees of freedom
    disps = np.zeros((tdof, 1))
    force = np.zeros((tdof, 1))
    sigma = np.zeros((numElem, 1))
    stiffness = np.zeros((tdof, tdof))
    np.set_printoptions(precision=3)

    # applying the load

    force[3] = -20000.0
    force[7] = -50000.0
    force[11] = -30000.0

    # defined boundary
    presDof = np.array([0, 1, 9])
    for e in range(numElem):
        indice = elemNodes[e, :]
        elemDof = np.array([indice[0] * 2, indice[0] * 2 + 1, indice[1] * 2,
                            indice[1] * 2 + 1])  # dof corresponding to the global stifnessmatrix
        xa = xx[indice[1]] - xx[indice[0]]  # length of the elemnts in x dir
        ya = yy[indice[1]] - yy[indice[0]]  # lenth of the elemnt in y dir
        len_elem = np.sqrt(xa * xa + ya * ya)
        c = xa / len_elem
        s = ya / len_elem
        k1 = (modE * Area_elem[e] / len_elem) * np.array([[c * c, c * s, -c * c, -c * s],
                                         [c * s, s * s, -c * s, -s * s],
                                         [-c * c, -c * s, c * c, c * s],
                                         [-c * s, -s * s, c * s, s * s]])
        stiffness[np.ix_(elemDof, elemDof)] += k1  # add elem K to the global K

    actDof = np.setdiff1d(np.arange(tdof), presDof)

# Correct way
    disp1 = np.linalg.solve(stiffness[np.ix_(actDof, actDof)], force[np.ix_(actDof)])


    disps[np.ix_(actDof)] = disp1.reshape([9,1])

    # stresses at elements

    for e in range(numElem):
        indice = elemNodes[e, :]
        elemDof = np.array([indice[0] * 2, indice[0] * 2 + 1, indice[1] * 2, indice[1] * 2 + 1])
        xa = xx[indice[1]] - xx[indice[0]]
        ya = yy[indice[1]] - yy[indice[0]]
        len_elem = np.sqrt(xa * xa + ya * ya)
        c = xa / len_elem
        s = ya / len_elem
        sigma[e] = (modE / len_elem) * np.dot(np.array([-c, -s, c, s]), disps[np.ix_(elemDof)])

    # print (disps)
    #print (sigma)
    # computing internal reactions
#    react = np.dot(stiffness, disps)
#    print (react.reshape((numNodes, 2)))
    for i,sig in enumerate(sigma):
        print("Elem: ",i, "\t Force: \t",sig*Area_elem[i],sig,sigma[i],Area_elem[i])
#    input("break")
    return np.max(np.abs(sigma))
def constraint1(A):
sum = 100000
for i in range(num_elem):
    sum = sum - A[i]
return sum


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

for i in range(200):
    print("Number iteration: ",i)
    sol = minimize(calculate_truss_force, x0, method='SLSQP', bounds=bnds,tol=0.0000000000001, constraints=con1)
    x0 = sol.x
    print(x0)
    print(calculate_truss_force(x0))

Answer 1

Gekko 使用自动微分为基于梯度的求解器提供导数。而不是在函数内部使用线性求解器:

disp1 = np.linalg.solve(stiffness[np.ix_(actDof, actDof)],\
                        force[np.ix_(actDof)])

这需要作为隐式 A x = b 而不是 x= A^-1 b 提供给 Gekko。

A = stiffness[np.ix_(actDof, actDof)]
b = np.array(force[np.ix_(actDof)])
disp1 = np.dot(A,b)

Gekko 同时求解方程，而不是像内部循环那样按顺序求解。此外，Gekko 仅对目标函数求值一次以构建符号表达式，然后将其编译为求解器的字节代码。它没有对 calculate_truss_force 的回调以对其进行多次评估。

另一个变化是使用 np.max 和 np.abs 到 m.max2 或 的 Gekko 版本m.max3 以及 m.min2 或 m.min3。这些是具有连续的一阶和二阶导数的 max 和 abs 的版本。 ...2 版本是 MPCC 而 ...3 版本是混合整数问题。

    for e in range(numElem):
        indice = elemNodes[e, :]
        elemDof = np.array([indice[0] * 2, indice[0] * 2 + 1, \
                            indice[1] * 2, indice[1] * 2 + 1])
        xa = xx[indice[1]] - xx[indice[0]]
        ya = yy[indice[1]] - yy[indice[0]]
        len_elem = np.sqrt(xa * xa + ya * ya)
        c = xa / len_elem
        s = ya / len_elem
        sigma[e] = m.abs2(m.Intermediate((modE / len_elem) * \
                          np.dot(np.array([-c, -s, c, s]), \
                          disps[np.ix_(elemDof)])))
    return m.max2(sigma)

如果您不想将目标函数用作“黑盒”，那么我推荐使用 scipy.optimize.minimize 而不是 Gekko。这是关于 Gekko and Scipy for optimization 的教程.