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

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

尝试:
从gekko进口GEKKO
除了:
＃pip安装gekko
进口点
pip.main(['install'，'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))

最佳答案

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来多次评估它。

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

    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)

如果您不想将目标函数用作“黑匣子”，那么我建议使用...3而不是Gekko。这是有关Gekko and Scipy for optimization的教程。

关于python-3.x - 尝试使用GEKKO OPTIMIZER时“No loop matching the specified signature and casting was found for ufunc solve”，我们在Stack Overflow上找到一个类似的问题：https://stackoverflow.com/questions/59865867/