尝试使用 GEKKO OPTIMIZER 时“找不到与指定签名和转换匹配的循环以进行 ufunc 解决”

Question

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

尝试： 从壁虎进口 GEKKO 除了： # pip install gekko 进口点子 pip.main(['安装','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 提供给Gekko，而不是x= A^-1 b。

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 用于Gekko 版本的m.max2 或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的教程。