gistfile1.txt

from math import pi
import numpy
from scipy.optimize import minimize
import context
from pystran import model
from pystran import section
from pystran import freedoms
from pystran import geometry
from pystran import plots

# The constraints are on the maximum deflection and the lowest frequency. The
# maximum deflection should be lower than the limit, and the fundamental
# frequencies should be above the lowest limit.
MAXTIPD = 10
LOWESTFREQ = 13

# The design variables are areas. A manufacturing constraint limits the
# smallest area.
INITIAL_AREA = pi * 60 * 7
SMALLEST_AREA = INITIAL_AREA / 10


# This function defines the model of the structure, based on the values of the
# design variables. The design variables are multipliers of the initial area,
# and so their values are around 1.0.
def tcant_model(dvs):
    m = model.create(2)
    model.add_joint(m, 1, [7000.0, 2500])
    model.add_joint(m, 2, [6000.0, 3500])
    model.add_joint(m, 3, [7000, -1000])
    model.add_joint(m, 4, [6000, 1500])
    model.add_joint(m, 5, [0, -1000])
    model.add_joint(m, 6, [0, 1500])
    model.add_support(m["joints"][5], freedoms.TRANSLATION_DOFS)
    model.add_support(m["joints"][6], freedoms.TRANSLATION_DOFS)
    W = 6000
    model.add_load(m["joints"][1], freedoms.U2, -W)
    model.add_load(m["joints"][2], freedoms.U2, -W)
    addM = 0.096
    model.add_mass(m["joints"][1], freedoms.U1, addM)
    model.add_mass(m["joints"][1], freedoms.U2, addM)
    E = 70000
    rho = 2.7e-9
    bar_connectivities = [
        [5, 3],
        [3, 1],
        [6, 4],
        [4, 2],
        [1, 2],
        [3, 4],
        [1, 4],
        [3, 6],
    ]
    for k, c in enumerate(bar_connectivities):
        s = section.truss_section(f"s{k}", E=E, A=INITIAL_AREA * dvs[k], rho=rho)
        model.add_truss_member(m, k, c, s)
    return m


# Helper function to compute the current mass of the structure.
def current_mass(m, dvs):
    mass = 0.0
    for member in m["truss_members"].values():
        A = INITIAL_AREA * dvs[member["mid"]]
        c = member["connectivity"]
        i, j = m["joints"][c[0]], m["joints"][c[1]]
        sect = member["section"]
        e_x, _, h = geometry.member_2d_geometry(i, j)
        mass += A * h * sect["rho"]
    return mass


# Helper function to compute the maximum deflection.
def max_tip_deflection(m):
    mtd = 0.0
    for j in m["joints"].values():
        mtd = max(mtd, max(abs(j["displacements"])))
    return mtd


# Helper function to compute the design responses. Both the statics and the
# free vibration problems need to be solved.
def solve(dvs):
    m = tcant_model(dvs)
    model.number_dofs(m)
    model.solve_statics(m)
    model.solve_free_vibration(m)
    drs = (current_mass(m, dvs), max_tip_deflection(m), m["frequencies"][0])
    return drs


# Initial guess -- these are the fractions of the initial area.
dvs0 = numpy.array([1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0])

# Report on the initial performance of the structure.
drs = solve(dvs0)
initial_mass = drs[0]
print("Initial areas: ", INITIAL_AREA * dvs0)
print("Initial mass: ", initial_mass)
print("Initial deflection: ", drs[1])
print("Initial frequency: ", drs[2])


# Objective function is the normalized mass.
def objective(dvs):
    drs = solve(dvs)
    return drs[0] / initial_mass


# Constraint function must return a non-negative value for each constraint.
# This is the constraint on the maximum deflection.
def constraint1(dvs):
    drs = solve(dvs)
    mtd = drs[1]
    return (MAXTIPD - mtd) / MAXTIPD


# Define a constraint on the frequency.
def constraint2(dvs):
    drs = solve(dvs)
    f = drs[2]
    return (f - LOWESTFREQ) / LOWESTFREQ


# Define constraints. They are both inequalities.
cons = [
    {"type": "ineq", "fun": constraint1},
    {"type": "ineq", "fun": constraint2},
]

# Define lower bounds for the design variables. There are no upper bounds.
bounds = [(SMALLEST_AREA / INITIAL_AREA, None) for _ in dvs0]

# Invoke the optimization function.
solution = minimize(
    objective,
    dvs0,
    method="SLSQP",
    bounds=bounds,
    constraints=cons,
    options={"ftol": 1e-7, "maxiter": 1000, "disp": True},
)

# Report on the solution of the optimization problem.
print(solution)

# Retrieve the values of the design variables, and compute the design responses
# for the optimal designed variables.
dvs = solution.x
drs = solve(dvs)

print("Optimal areas: ", INITIAL_AREA * dvs)
print("Optimal mass: ", drs[0])
print("Optimal deflection: ", drs[1])
print("Optimal frequency: ", drs[2])