Scipy example

[germa89]

Take this example and make it part of the documentation, ideally comparing its results with a trust made with beam elements and/or link elements.

Truss Example

Problem Description:

We are tasked with optimizing the cross-sectional areas of the members of a truss structure to minimize its total weight, while ensuring that the stress in each member does not exceed a predefined limit (i.e., yield stress). This involves solving a constrained optimization problem, where the constraints are based on structural safety criteria.

Model:

A truss is a structure made up of straight members connected at joints (nodes).
We minimize the total weight of the truss by adjusting the cross-sectional areas of the members.
The constraints ensure that stresses in the members due to applied loads are within safe limits.
The stress in each member is computed using the force-displacement relationship based on Hooke’s Law.

Assumptions:

• We will use linear elasticity.
• Small deformations, so linear analysis applies.
• The truss is statically determinate (i.e., forces in members are calculated using equilibrium equations).

Steps:

1. Define the truss geometry (nodes and members).
2. Define the force-displacement equations.
3. Set the objective function: minimize total weight.
4. Set the stress constraints for each member.
5. Use SciPy’s optimization tools to solve the problem.

Code

import numpy as np
from scipy.optimize import minimize


###############################################################################
## Problem definition
max_displacement = 0.001
minimal_section = 0.001

# Density of the material (steel)
density = 7850  # kg/m^3

# Young's Modulus (material property)
E = 210e9  # Pa (steel)

# Stress constraint: stress in each member must be below yield stress
yield_stress = 250e6  # Pa (steel yield stress)

# Geometry of the truss bridge (nodes and members)
nodes = np.array([
    [0, 0],  # Node 0: Pinned support
    [2, 0],  # Node 1
    [4, 0],  # Node 2 - Load -15 vertical
    [6, 0],  # Node 3
    [8, 0],  # Node 4 - Load -15 vertical
    [10, 0], # Node 5: Roller support
    [2, 3],  # Node 6
    [4, 3],  # Node 7
    [6, 3],   # Node 8
    [8, 3]   # Node 9
])

members = [
    (0, 1), (1, 2), (2, 3), (3, 4), (4, 5), # Lower members
    (0, 6), (1, 7), (2, 8), (3, 7), (4, 8), (5, 9), # Diagonal members
    (1, 6), (2, 7), (3, 8), (4, 9), # vertical members
    (6, 7), (7, 8), (8, 9), # Upper members
]

# Number of members and nodes
n_members = len(members)
n_nodes = len(nodes)

# Applied loads at each node (Node 2 and Node 3 experience loads)
loads = np.zeros((n_nodes, 2))  # No load on other nodes
loads[2] = [0, -15]  # Vertical load (-150 kN) on Node 2
loads[3] = [10, -15]  # Horizontal (5 kN) and vertical (-150 kN) loads on Node 3

multiplier=1E4 # KiloNewtons
loads = loads*multiplier

# Displacements
fixed_dofs = [0, 1, 11]  # Node 0 and 5 fixed in both directions

###############################################################################
## Plotting function

import matplotlib.pyplot as plt
import numpy as np

def plot_truss(nodes, members, areas=None, loads=None, constraints=None, displacements=None, title="Truss Structure"):
    """
    Plots the truss structure using node coordinates and member connections.
    
    Parameters:
        nodes (np.array): An array of shape (n_nodes, 2) representing the (x, y) coordinates of the nodes.
        members (list of tuples): A list of tuples where each tuple contains two node indices representing a member.
        areas (np.array or None): An array of cross-sectional areas for each member. Default is None (uniform thickness).
        title (str): The title of the plot.
    """
    plt.figure(figsize=(8, 6))
    font = {'weight': 'bold'}
    plt.title(title, fontdict=font)
    
    # Default line width if no areas provided
    if areas is None:
        areas = np.ones(len(members))
    
    # Normalize areas for line width visualization (for better plotting)
    min_thickness = 1
    max_thickness = 8

    if areas is not None and (areas.max() - areas.min()) != 0:
        norm_areas = (areas - areas.min()) / (areas.max() - areas.min())
        line_widths = norm_areas * (max_thickness - min_thickness) + min_thickness
    else:
        line_widths = np.ones(len(members)) * 2  # Default line width

    # Plot all the members
    for i, (n1, n2) in enumerate(members):
        x_coords = [nodes[n1][0], nodes[n2][0]]
        y_coords = [nodes[n1][1], nodes[n2][1]]
        plt.plot(x_coords, y_coords, 'b-', lw=line_widths[i])
    
    # Plot the nodes as red circles
    plt.scatter(nodes[:, 0], nodes[:, 1], color='r', zorder=5)
    
    # Label the nodes
    for i, (x, y) in enumerate(nodes):
        plt.text(x, y, f'{i}', fontsize=12, ha='right')
    

    # Plot loads as arrows
    if loads is not None:
        for i, (fx, fy) in enumerate(loads):
            if fx != 0 or fy != 0:
                plt.arrow(nodes[i][0], nodes[i][1], fx * 0.1/multiplier, fy * 0.1/multiplier, head_width=0.2, head_length=0.3, fc='g', ec='g')
    
    # Plot constraints (supports)
    if constraints is not None:
        for node in constraints:
            marker = 6 if node % 2 else 4
            node = node//2
            plt.scatter(nodes[node][0], nodes[node][1], marker=marker, color="black", s=100, zorder=10)

    # Set axis limits and labels
    plt.xlabel('X-coordinate (m)')
    plt.ylabel('Y-coordinate (m)')
    plt.gca().set_aspect('equal', adjustable='box')
    plt.grid(True)
    plt.show()

plot_truss(nodes, members, loads=loads, constraints=fixed_dofs, title="Truss Bridge Structure")

###############################################################################
## Solving
# Length and initial cross-sectional areas of the members
def member_length(n1, n2):
    return np.linalg.norm(nodes[n1] - nodes[n2])

lengths = np.array([member_length(n1, n2) for n1, n2 in members])

# Objective function: minimize weight of the truss
def weight(areas):
    w = np.sum(areas * lengths * density)
    # print(w)
    return w

def weight_opt(areas):
    # Normalize
    areas = areas/np.abs(areas)
    return weight(areas)

# Force-displacement matrix (stiffness matrix) and stress calculations
def stiffness_matrix(areas):
    K = np.zeros((2 * n_nodes, 2 * n_nodes))  # Global stiffness matrix
    for i, (n1, n2) in enumerate(members):
        L = lengths[i]
        A = areas[i]
        cos_theta = (nodes[n2][0] - nodes[n1][0]) / L
        sin_theta = (nodes[n2][1] - nodes[n1][1]) / L
        k_local = E * A / L * np.array([
            [ cos_theta**2,  cos_theta*sin_theta, -cos_theta**2, -cos_theta*sin_theta],
            [ cos_theta*sin_theta, sin_theta**2, -cos_theta*sin_theta, -sin_theta**2],
            [-cos_theta**2, -cos_theta*sin_theta,  cos_theta**2,  cos_theta*sin_theta],
            [-cos_theta*sin_theta, -sin_theta**2,  cos_theta*sin_theta, sin_theta**2]
        ])
        # Assemble the local stiffness matrix into the global stiffness matrix
        dof_map = [2*n1, 2*n1+1, 2*n2, 2*n2+1]
        for row in range(4):
            for col in range(4):
                K[dof_map[row], dof_map[col]] += k_local[row, col]
    return K

# Displacement constraint function (nodal equilibrium equations)
def displacements(areas):
    K = stiffness_matrix(areas)
    # Boundary conditions
    free_dofs = np.setdiff1d(np.arange(2 * n_nodes), fixed_dofs)

    # Solve for displacements at free nodes
    K_reduced = K[np.ix_(free_dofs, free_dofs)]
    load_reduced = loads.flatten()[free_dofs]
    u_free = np.linalg.solve(K_reduced, load_reduced)

    displacements_full = np.zeros(2 * n_nodes)
    displacements_full[free_dofs] = u_free
    return displacements_full.reshape((n_nodes, 2))

# Stress constraint
def stress_constraint(areas):
    # https://people.duke.edu/~hpgavin/cee421/truss-method.pdf
    u = displacements(areas)
    stresses = np.zeros(n_members)
    for i, (n1, n2) in enumerate(members):
        L = lengths[i]
        A = areas[i]
        cos_theta = (nodes[n2][0] - nodes[n1][0]) / L
        sin_theta = (nodes[n2][1] - nodes[n1][1]) / L

        delta = u[n2] - u[n1]  # Displacement difference between nodes
        force = E * A / L * (delta[0]*cos_theta + delta[1]*sin_theta)
        stresses[i] = abs(force / A)

    return yield_stress - stresses  # Positive means constraint satisfied

# Add displacement constraint
def displacement_constraint(areas):
    u = displacements(areas)
    return  max_displacement - (u[:,0]**2 + u[:,1]**2)**0.5

# Add positivity constraint for cross-sectional areas (they must be positive)
def area_positivity_constraint(areas):
    return (areas - minimal_section)   # Positive areas, minimal area

# Initial guess for cross-sectional areas (0.01 m^2 for each member)
initial_areas = np.ones(n_members) * minimal_section * 2

# Solve the optimization problem with constraints
constraints = [
    {'type': 'ineq', 'fun': stress_constraint},
    {'type': 'ineq', 'fun': displacement_constraint},
    {'type': 'ineq', 'fun': area_positivity_constraint}
]

# Solve
result = minimize(
        weight_opt,
        initial_areas,
        constraints=constraints,
        method="SLSQP", 
        options={"disp": True}
        )

# Optimized areas and corresponding weight
optimized_areas = result.x
optimized_weight = weight(optimized_areas)

# Results
print("Optimized Cross-Sectional Areas (m^2):", optimized_areas)
print("Optimized Truss Weight (kg):", optimized_weight)

# Plotting result
plot_truss(nodes, members, areas=optimized_areas, title="Result Truss Bridge Structure")

Optimization terminated successfully    (Exit mode 0)
            Current function value: 389621.46507435385
            Iterations: 7
            Function evaluations: 134
            Gradient evaluations: 7
Optimized Cross-Sectional Areas (m^2): [0.01898861 0.02461479 0.02264893 0.01983588 0.01218406 0.02175955
 0.02175955 0.00813483 0.01038501 0.02239731 0.02239731 0.01817839
 0.01332413 0.0154459  0.01875639 0.01176874 0.02131085 0.01218406]
Optimized Truss Weight (kg): 6824.634294298257