I am trying to optimize the weight of a truss structure by changing the areas. I have this (shortened/verified code for FEA truss) where x is a array of element areas:
stresses = np.zeros(10)
def FEA(x):
#code to calculate stresses, weight
stresses = element_stresses
return weight
I then have a constraint function like this:
def eval_constraints(x):
global stresses
constraints = [stresses[0]/25-1, stresses[0]/-25-1,
stresses[1]/25-1, stresses[1]/-25-1,
stresses[2]/25-1, stresses[2]/-25-1,
stresses[3]/25-1, stresses[3]/-25-1,
stresses[4]/25-1, stresses[4]/-25-1,
stresses[5]/25-1, stresses[5]/-25-1,
stresses[6]/25-1, stresses[6]/-25-1,
stresses[7]/25-1, stresses[7]/-25-1,
stresses[8]/75-1, stresses[8]/-75-1,
stresses[9]/25-1, stresses[9]/-25-1]
return constraints
The optimize command is as such:
x0 = np.full(10,5.)
bnds = [(0.0, None)] * 10
cons = [{'type': 'ineq', 'fun': eval_constraints}]
result = minimize(FEA, x0, bounds=bnds, constraints=cons)
What happens now is that the derivatives with respect to each x variable get computed. But then the optimzation suggets 0.0 for all the areas. This violates the constraints:
[5. 5. 5. 5. 5. 5. 5. 5. 5. 5.]
[5.00000001 5. 5. 5. 5. 5. 5. 5. 5. 5.]
[5. 5.00000001 5. 5. 5. 5. 5. 5. 5. 5.]
[5. 5. 5.00000001 5. 5. 5. 5. 5. 5. 5.]
[5. 5. 5. 5.00000001 5. 5. 5. 5. 5. 5.]
[5. 5. 5. 5. 5.00000001 5. 5. 5. 5. 5.]
[5. 5. 5. 5. 5. 5.00000001 5. 5. 5. 5.]
[5. 5. 5. 5. 5. 5. 5.00000001 5. 5. 5.]
[5. 5. 5. 5. 5. 5. 5. 5.00000001 5. 5.]
[5. 5. 5. 5. 5. 5. 5. 5. 5.00000001 5.]
[5. 5. 5. 5. 5. 5. 5. 5. 5. 5.00000001]
[0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]
I even try setting bnds = [(0.1, None)] * 10, but it doesn't move that far from the x0 value.
'Positive directional derivative for linesearch' error
[4.99996444, 4.99996444, 4.99996444, 4.99996444, 4.99996444,
4.99996444, 4.99998758, 4.99997551, 4.99996445, 5.00053543]
The optimized value is [7.90, 0.10, 8.10, 3.90, 0.10, 0.10, 5.80, 5.51, 3.68, 0.14]. Any help would be appreciated. Thanks in advance.
import numpy as np
from numpy import cos,sin
from scipy.optimize import minimize, fmin_slsqp
#nodes: x, y
nodes = np.array([[720, 360],
[720, 0],
[360,360],
[360,0],
[0,360],
[0,0]])
#materials: E, rho, sigmaT, sigmaC
materials = np.array([[10000, 0.1, 25, 25],
[10000, 0.1, 75, 75]])
#boundary conditions: nodeID, xDOF, yDOF
boundaries = np.array([[5, 0, 0],
[6, 0, 0]])
#concentrated loads: nodeID, xLoad, yLoad
cloads = np.array([[2, 0, -100],
[4, 0, -100]])
#generic element stiffness matrix
bar_truss_element_k = np.array([[1, 0, -1, 0],
[0, 0, 0, 0],
[-1, 0, 1, 0],
[0, 0, 0, 0]])
def T_matrix(theta):
T = np.zeros((4,4))
T[0,0] = T[1,1] = T[2,2] = T[3,3] = cos(theta)
T[0,1] = T[2,3] = sin(theta)
T[1,0] = T[3,2] = -sin(theta)
return T
def element_L_and_theta(node1, node2):
diff = node1-node2
L = np.linalg.norm(diff)
theta = np.arctan2(diff[1],diff[0])
return L, theta
def FEA(x):
#elements: node1, node2, A, matID
elements = np.array([[3, 5, x[0], 1],
[1, 3, x[1], 1],
[4, 6, x[2], 1],
[2, 4, x[3], 1],
[3, 4, x[4], 1],
[1, 2, x[5], 1],
[4, 5, x[6], 1],
[3, 6, x[7], 1],
[2, 3, x[8], 2],
[1, 4, x[9], 1]])
#Intialize matricies and vectors
num_DOFs = nodes.shape[0] * 2
stiffness_matrix = np.zeros((num_DOFs, num_DOFs))
forces = np.zeros((num_DOFs,1))
for nodeID, xLoad, yLoad in cloads:
nodeID = int(nodeID)
forces[2*(nodeID-1)] = xLoad if xLoad != 0 else 0
forces[2*(nodeID-1)+1] = yLoad if yLoad != 0 else 0
#Compute the stiffness matrix
for element in elements:
node1, node2, A, matID = element
node1, node2, matID = int(node1)-1, int(node2)-1,
int(matID)-1
E = materials[matID][0]
L, theta = element_L_and_theta(nodes[node1],nodes[node2])
T = T_matrix(theta)
element_k = (A*E/L) * np.transpose(T) @
bar_truss_element_k @ T
stiffness_matrix[np.ix_([2*node1, 2*node1+1, 2*node2,
2*node2+1], [2*node1, 2*node1+1, 2*node2, 2*node2+1])] +=
element_k
#Identify and eliminate zero DOFs
zero_DOFs = [2*(condition[0]-1)+i for condition in boundaries
for i,DOF in enumerate(condition[1:]) if DOF==0]
non_zero_mask = np.ones(num_DOFs, dtype=bool)
non_zero_mask[zero_DOFs] = False
reduced_stiffness_matrix = stiffness_matrix[non_zero_mask][:,
non_zero_mask]
reduced_forces = forces[non_zero_mask]
#Solve for the displacements
reduced_displacements =
np.linalg.solve(reduced_stiffness_matrix, reduced_forces)
displacements = np.zeros((num_DOFs,1))
displacements[np.setdiff1d(np.arange(num_DOFs), zero_DOFs)] =
reduced_displacements
#Solve for the forces
forces = stiffness_matrix @ displacements
#Solve for the element stresses
element_stresses = np.zeros(elements.shape[0])
weight = 0.0
for i, element in enumerate(elements):
node1, node2, A, matID = element
node1, node2, matID = int(node1)-1, int(node2)-1, int(matID)-1
E = materials[matID][0]
L, theta = element_L_and_theta(nodes[node1],nodes[node2])
element_stresses[i] = (E / L) * np.array([np.cos(theta),
np.sin(theta), -np.cos(theta), -np.sin(theta)]) @
displacements[[2*node1, 2*node1+1, 2*node2, 2*node2+1]]
weight += L*A*0.1 #rho=0.1
stresses = element_stresses
return weight
def eval_constraints(x):
global stresses
constraints = [1-stresses[0]/25, 1-stresses[0]/-25,
1-stresses[1]/25, 1-stresses[1]/-25,
1-stresses[2]/25, 1-stresses[2]/-25,
1-stresses[3]/25, 1-stresses[3]/-25,
1-stresses[4]/25, 1-stresses[4]/-25,
1-stresses[5]/25, 1-stresses[5]/-25,
1-stresses[6]/25, 1-stresses[6]/-25,
1-stresses[7]/25, 1-stresses[7]/-25,
1-stresses[8]/75, 1-stresses[8]/-75,
1-stresses[9]/25, 1-stresses[9]/-25]
return constraints
x0 = np.full(10,5.)
FEA(x0)
#x_correct = np.array([7.90, 0.10, 8.10, 3.90, 0.10, 0.10, 5.80,
5.51, 3.68, 0.14])
bnds = [(0.1, None)] * 10
cons = [{'type': 'ineq', 'fun': eval_constraints}]
result = minimize(FEA, x0, bounds=bnds, constraints=cons)
print(result)