«Не найдена l oop, соответствующая указанной подписи, и преобразование было найдено для ufun c solve» при попытке использовать GEKKO OPTIMIZER - PullRequest
Новая
       6

«Не найдена l oop, соответствующая указанной подписи, и преобразование было найдено для ufun c solve» при попытке использовать GEKKO OPTIMIZER

3 голосов
/ 22 января 2020

Я пытаюсь оптимизировать свое решение для фермы FEA, используя GEKKO. Поэтому я определил все решения FEA как одну функцию и попытался сделать силы переменными. Так что цель состоит в том, чтобы минимизировать силы одним ограничением (сумма всех сил = -100000). Но после запуска кода я получаю ошибку, которая Я не мог исправить это, поиск в Google также не помог (я новичок в python). Буду признателен, если кто-то может сказать мне, где проблема. Заранее спасибо.

Попробовать: из gekko импортировать GEKKO кроме: # pip install gekko import pip pip.main (['install', 'gekko']) из gekko import 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))

1 Ответ

2 голосов
/ 23 января 2020

Gekko использует автоматическое дифференцирование c для предоставления производных для решателей на основе градиента. Вместо использования линейного решателя внутри вашей функции: 

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

это нужно передать Gekko как неявное A x = b, а не x= A^-1 b. 

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

Gekko одновременно решает уравнение, а не последовательно, как с внутренним l oop. Кроме того, Gekko оценивает целевую функцию только один раз, чтобы построить выражения Symboli c, которые затем компилируется в байт-код для решателя. У него нет обратного вызова на calculate_truss_force для его оценки несколько раз.

Еще одно изменение - использование np.max и np.abs для версий Gekko m.max2 или m.max3, а также с m.min2 или m.min3. Это версии max и abs, которые имеют непрерывные первую и вторую производные. Версия ...2 является MP CC, а версия ...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 и Scipy для оптимизации .

...