Tiawichi Research

Structural optimization has been an active field of research for many years. With the advent of finite element analysis methods, the promise of structural optimization became even brighter, especially in weight sensitive structures such as aerospace vehicles. The initial focus of the methods was on gradient methods in which the derivatives of the response WRT the variables being optimized was required (Vanderplaats). Within the last decade or two, less demanding methods have been proposed and coupled with the wide availability of compute clusters, methods for effective structural optimization methodologies are realizable.

Some of the tools are called by various names such as the method of simulated annealing, ant colony, genetic algorithms, parasitic swarm optimization, evolutionary structural optimization. Several of these methods as well as more classical methods have been implemented in the scipy library.

In the course of delivering FEA solutions to our clients, we often run into constrained optimization problems that are not very tractable with the traditional gradient methods. A quick solution was to couple a scipy method with the non-linear finite element analysis as such:

import os,sys
import glob
from scipy.optimize import cobyla as co
import SEL

job_id = "Dam_shell_2"

[FLANGE,
TABS,
TABS_THIN,
RIM_FLANGE,
CIRCUM,
PLATE,
WEB,
RIM,
RIM_THICK,
FLANGE_RADIAL]=range(10)

# Initial guess for solution
x_init=[]
x_init.insert(FLANGE,0.3)
x_init.insert(TABS,1.0)
x_init.insert(TABS_THIN,0.4)
x_init.insert(RIM_FLANGE,0.3125)
x_init.insert(CIRCUM,0.31)
x_init.insert(PLATE,0.165)
x_init.insert(WEB,0.165)
x_init.insert(RIM,0.33)
x_init.insert(RIM_THICK,0.86)
x_init.insert(FLANGE_RADIAL,0.15)

def weight(x):
   '''Objective function to minimize'''

    flange=max(0.001,x[FLANGE])
    tabs=max(0.001,x[TABS])
    tabs_thin=max(0.001,x[TABS_THIN])
    rim_flange=max(0.001,x[RIM_FLANGE])
    circum=max(0.001,x[CIRCUM])
    plate=max(0.001,x[PLATE])
    web=max(0.001,x[WEB])
    rim=max(0.001,x[RIM])
    rim_thick=max(0.001,x[RIM_THICK])
    flange_radial=max(0.001,x[FLANGE_RADIAL])

    # Modify parameters using the SEL package to replace the parameters in the
    # FEA input file with the modified parameters

    params = '~flange=[0123456789.]*~=flange\=%f '%flange
    params += '~tabs=[0123456789.]*~=tabs\=%f '%tabs
    params += '~tabs_thin=[0123456789.]*~=tabs_thin\=%f '%tabs_thin
    params += '~rim_flange=[0123456789.]*~=rim_flange\=%f '%rim_flange
    params += '~circum=[0123456789.]*~=circum\=%f '%circum
    params += '~plate=[0123456789.]*~=plate\=%f '%plate
    params += '~web=[0123456789.]*~=web\=%f '%web
    params += '~rim=[0123456789.]*~=rim\=%f '%rim
    params += '~rim_thick=[0123456789.]*~=rim_thick\=%f '%rim_thick
    params += '~flange_radial=[0123456789.]*~=flange_radial\=%f '%flange_radial

    parameters=SEL.SEL(params)
    parameters.set(file_handling_flag = 3)
    parameters(job_id+'.inp')

    # Clean up files from previous analyses
    files = glob.glob(job_id+".*")
    for file in files:
        if file[-3:] != "inp":
            os.remove(file)

    # Execute FEA analysis
    os.spawnv(os.P_WAIT,"<path to fea code>",["args"])

    # Extract and return weight
    os.spawnv(os.P_WAIT,"<path to weight extractor>",["args"])

    f = open(job_id+"_Weight.txt")
    w = float(f.readline())
    return w

def stress(x):
    '''Constraint function on stress'''

    # Extract and return maximum stress
    os.spawnv(os.P_WAIT,"<path to stress extractor>",["args"])
    f = open(job_id+"_MaxMises.txt")
    sigma = float(f.readline())
    # Normalize constraint, 28000.0 is the allowable stress
    return 28000.0/sigma-1.0

if __name__ == '__main__':
    '''Minimize a function using the Constrained Optimization BY Linear
    Approximation (COBYLA) method'''
    print co.fmin_cobyla(weight,
                         x_init,
                         [stress],
                         rhobeg=0.01,
                         rhoend=0.0001,
                         iprint=3,
                         maxfun=300)