Tutorial: discrete-based lightweight optimization using the ghost layer method#
Date: 31 of August 2024
Author: Saullo G. P. Castro
Cite this tutorial as:
Castro, SGP. Methods for analysis and design of composites (Version 0.7.1) [Computer software]. 2024. https://doi.org/10.5281/zenodo.2871782
Installing required modules#
[115]:
!python -m pip install numpy scipy pymoo composites > tmp.txt
DEPRECATION: Loading egg at /Users/saullogiovanip/miniconda3/lib/python3.12/site-packages/composites-0.7.0-py3.12-macosx-11.1-arm64.egg is deprecated. pip 24.3 will enforce this behaviour change. A possible replacement is to use pip for package installation.. Discussion can be found at https://github.com/pypa/pip/issues/12330
DEPRECATION: Loading egg at /Users/saullogiovanip/miniconda3/lib/python3.12/site-packages/fea-0.0.1-py3.12.egg is deprecated. pip 24.3 will enforce this behaviour change. A possible replacement is to use pip for package installation.. Discussion can be found at https://github.com/pypa/pip/issues/12330
DEPRECATION: Loading egg at /Users/saullogiovanip/miniconda3/lib/python3.12/site-packages/docopt-0.6.2-py3.12.egg is deprecated. pip 24.3 will enforce this behaviour change. A possible replacement is to use pip for package installation.. Discussion can be found at https://github.com/pypa/pip/issues/12330
DEPRECATION: Loading egg at /Users/saullogiovanip/miniconda3/lib/python3.12/site-packages/coverage-7.5.4-py3.12-macosx-11.1-arm64.egg is deprecated. pip 24.3 will enforce this behaviour change. A possible replacement is to use pip for package installation.. Discussion can be found at https://github.com/pypa/pip/issues/12330
DEPRECATION: Loading egg at /Users/saullogiovanip/miniconda3/lib/python3.12/site-packages/coverage-7.6.0-py3.12-macosx-11.1-arm64.egg is deprecated. pip 24.3 will enforce this behaviour change. A possible replacement is to use pip for package installation.. Discussion can be found at https://github.com/pypa/pip/issues/12330
DEPRECATION: Loading egg at /Users/saullogiovanip/miniconda3/lib/python3.12/site-packages/setuptools_git_version-1.0.3-py3.12.egg is deprecated. pip 24.3 will enforce this behaviour change. A possible replacement is to use pip for package installation.. Discussion can be found at https://github.com/pypa/pip/issues/12330
DEPRECATION: Loading egg at /Users/saullogiovanip/miniconda3/lib/python3.12/site-packages/structsolve-0.2.2-py3.12.egg is deprecated. pip 24.3 will enforce this behaviour change. A possible replacement is to use pip for package installation.. Discussion can be found at https://github.com/pypa/pip/issues/12330
DEPRECATION: Loading egg at /Users/saullogiovanip/miniconda3/lib/python3.12/site-packages/coveralls-4.0.1-py3.12.egg is deprecated. pip 24.3 will enforce this behaviour change. A possible replacement is to use pip for package installation.. Discussion can be found at https://github.com/pypa/pip/issues/12330
DEPRECATION: Loading egg at /Users/saullogiovanip/miniconda3/lib/python3.12/site-packages/pyfe3d-0.5.0-py3.12-macosx-11.1-arm64.egg is deprecated. pip 24.3 will enforce this behaviour change. A possible replacement is to use pip for package installation.. Discussion can be found at https://github.com/pypa/pip/issues/12330
DEPRECATION: Loading egg at /Users/saullogiovanip/miniconda3/lib/python3.12/site-packages/tudaesasII-2024.3-py3.12.egg is deprecated. pip 24.3 will enforce this behaviour change. A possible replacement is to use pip for package installation.. Discussion can be found at https://github.com/pypa/pip/issues/12330
References#
Riche, R. Le, and Haftka, R. T., 1993, “Optimization of Laminate Stacking Sequence for Buckling Load Maximization by Genetic Algorithm,” AIAA J., 31(5), pp. 951–956. https://arc.aiaa.org/doi/10.2514/3.11710
Castro, S. G. P., and Guimar, T. A. M., 2010, “Comparison of Free Stacking Sequence Approach versus a Predefined 0 / + 45 / -45 / 90 Sequence in a Typical Aircraft Wing Optimization,” 2nd International Conference on Engineering Optimization (ENGOPT), Lisbon, Portugal, pp. 1–11. https://www.researchgate.net/publication/233406751_Comparison_of_free_stacking_sequence_approach_versus_a_predefined_045-4590_sequence_in_a_typical_aircraft_wing_optimization
Schläpfer, B., and Kress, G., 2012, “A Sensitivity-Based Parameterization Concept for the Automated Design and Placement of Reinforcement Doublers,” Compos. Struct., 94(3), pp. 896–903. https://linkinghub.elsevier.com/retrieve/pii/S0263822311003655
Ghost layer, or free stacking sequence, or “T-THETA”, …
based on the same principle of creating dummy layers that can removed or added during the optimization
Main idea behind the ghost layer method:#
duplicate the number of design variables to be able to add or remove layers
Laminate plate geometry and applied loading:#
Constraint - Critical buckling load \(\lambda_b\), calculated analytically with:#
From Riche and Haftka’s paper:
\[\frac{\lambda_b(m, n)}{\pi^2} = \frac{\left[D_{11} (m/a)^4 + 2 (D_{12} + 2 D_{66}) (m/a)^2 (n/b)^2 + D_{22} (n/b)^4\right]}{((m/a)^2 Nxx + (n/b)^2 Nyy)}\ (1)\]
Constraint - strain failure#
This constraint requires that all strain components remain below an allowable limit, assuming that \(\gamma_{xy}\) will be zero for this bi-axial loading, given that the laminate is balanced.
The principal strains on each layer can be calculated solving first the laminate bi-axial strains \(\epsilon_x\) and \(\epsilon_y\); and later the in-plane strains of each layer: \(\epsilon_1^i\), \(\epsilon_2^i\) and \(\gamma_{12}^i\):
[2]:
# NOTE allowable strains from the reference paper
epsilon_1_allowable = 0.008
epsilon_2_allowable = 0.029
gamma_12_allowable = 0.015
Optimization description#
Take as reference the first design case of Table 3
Note that the design is currently constrained by the failure index.
Note that the number of plies here has been pre-defined to 48 plies, and with the ghost layer approach, we will be able to remove plies.
We will set a target design load by means of a target lambda that will multiply the bi-axial loading state.
Defining the reference case#
Using the first stacking sequence of the Table 3 above:
[3]:
# Table 03, first row
a_value = 20. # [in]
b_value = 5. # [in]
load_Nxx_unit = 1. # [lb]
load_Nyy_unit = 0.125 # [lb]
stack_ref_half = [90]*2 + [+45,-45]*4 + [0]*4 + [+45, -45] + [0]*4 + [+45, -45] + [0]*2
stack_ref = stack_ref_half + stack_ref_half[::-1] # symmetry condition
assert len(stack_ref) == 48
num_variables = len(stack_ref) // 4 # NOTE independent angles for a symmetric and balanced laminate
# NOTE make sure that you understand this
num_variables = 2*num_variables # NOTE twice as many variables needed for the ghost layer approach
# material properties
E11 = 127.55e9 # [Pa]
E22 = 13.03e9
nu12 = 0.3
G12 = 6.41e9
G13 = 6.41e9
G23 = 6.41e9
laminaprop = (E11, E22, nu12, G12, G13, G23)
ply_thickness = 0.005*25.4/1000 # [m]
h = ply_thickness*len(stack_ref)
[4]:
from composites import laminated_plate
prop_ref = laminated_plate(stack=stack_ref, plyt=ply_thickness, laminaprop=laminaprop)
Analytical equation to calculate buckling#
[5]:
from numpy import pi
def calc_buckling_analytical(prop):
a = a_value*25.4/1000 # [m] along x
b = b_value*25.4/1000 # [m] along y
Nxx = load_Nxx_unit*4.448222/(25.4/1000) #[N/m]
Nyy = load_Nyy_unit*4.448222/(25.4/1000) #[N/m]
D11 = prop.ABD[3, 3]
D12 = prop.ABD[3, 4]
D22 = prop.ABD[4, 4]
D66 = prop.ABD[5, 5]
lambda_b_min = 1e30
for m in range(1, 21):
for n in range(1, 21):
lambda_b = (pi**2*(D11*(m/a)**4
+ 2*(D12 + 2*D66)*(m/a)**2*(n/b)**2
+ D22*(n/b)**4)/((m/a)**2*Nxx + (n/b)**2*Nyy)
)
lambda_b_min = min(lambda_b_min, lambda_b)
return lambda_b_min
Function to calculate strain failure criterion#
[6]:
import numpy as np
import scipy.optimize as opt
def calc_failure_load_Haftka(prop):
def strain_MS(lbd):
Nxx = -load_Nxx_unit*4.448222/(25.4/1000) #[N/m]
Nyy = -load_Nyy_unit*4.448222/(25.4/1000) #[N/m]
Nxx = lbd*Nxx*1.5
Nyy = lbd*Nyy*1.5
vecN = np.asarray([Nxx, Nyy])
A11 = prop.ABD[0, 0]
A12 = prop.ABD[0, 1]
A22 = prop.ABD[1, 1]
exx, eyy = np.linalg.inv(np.array([[A11, A12], [A12, A22]])) @ vecN
margin_of_safety = 1e15
for thetadeg in prop.stack:
cost = np.cos(np.deg2rad(thetadeg))
sint = np.sin(np.deg2rad(thetadeg))
epsilon_i_1 = cost**2*exx + sint**2*eyy
epsilon_i_2 = sint**2*exx + cost**2*eyy
gamma_i_12 = sint**2*(eyy - exx)
ms_new = min(
epsilon_1_allowable/abs(epsilon_i_1) - 1,
epsilon_2_allowable/abs(epsilon_i_2) - 1,
gamma_12_allowable/abs(gamma_i_12) - 1
)
margin_of_safety = min(margin_of_safety, ms_new)
return margin_of_safety
lbd_init = 100
positiveMS = opt.NonlinearConstraint(strain_MS, 0., np.inf, jac='2-point')
res = opt.minimize(strain_MS, lbd_init, tol=1e-6, bounds=((100, None),),
constraints=[positiveMS], jac='2-point')
assert res.success
return res.x[0]
Verifying constraint functions#
Compare with columns 2 and 3 of Table 3 above, first row.
[7]:
lambda_cs_ref = calc_failure_load_Haftka(prop_ref)
print('lambda_cs_ref', lambda_cs_ref)
assert np.isclose(lambda_cs_ref, 13518.66, rtol=1e-3)
lambda_cb_ref_analytical = calc_buckling_analytical(prop_ref)
print('lambda_cb_ref_analytical', lambda_cb_ref_analytical)
lambda_cs_ref 13517.970921520578
lambda_cb_ref_analytical 14167.523183899806
/var/folders/x7/jww7lc053bv10s9m981_324dmy7lhy/T/ipykernel_57190/2488278243.py:25: RuntimeWarning: divide by zero encountered in divide
gamma_12_allowable/abs(gamma_i_12) - 1
Defining design load#
[79]:
target_lambda = 10000
Stacking sequence and thickness from integer vector#
Half of the \(x\) vector constains information about the layup, the other half about the corresponding ply thickness, implementing the ghost layer approach.
[80]:
available_angles = [0, 45, 90]
ply_thickness = 0.005*25.4/1000 # [m]
def layup_thick_from_x(x):
layup = []
plyts = []
for xi in x[0:len(x)//2]:
layup.append(available_angles[xi])
for xi in x[len(x)//2:]:
plyts.append(ply_thickness*xi)
layup_bal = []
plyts_bal = []
for angle_deg, plyt in zip(layup, plyts): # balancing
layup_bal.append(angle_deg)
plyts_bal.append(plyt)
if angle_deg != 0 and angle_deg != 90:
layup_bal.append(-angle_deg)
plyts_bal.append(plyt)
else:
layup_bal.append(angle_deg)
plyts_bal.append(plyt)
layup_bal = layup_bal + layup_bal[::-1] # symmetry
plyts_bal = plyts_bal + plyts_bal[::-1] # symmetry
non_zero = np.logical_not(np.asarray(plyts_bal) == 0)
return np.asarray(layup_bal)[non_zero], np.asarray(plyts_bal)[non_zero]
Defining objective function#
[81]:
#def objective_Workshop1(x): # NOTE to be minimized
# stack = discrete_stack_from_continuous_x(x)
# prop = laminated_plate(stack=stack, plyt=ply_thickness, laminaprop=laminaprop)
# lambda_cb = calc_buckling_FE(prop)
# lambda_cs = calc_failure_load_Haftka(prop)
# p = 0.08 # NOTE from the reference paper
# obj = (1 - p)*min(lambda_cs, lambda_cb)
# return 1/obj
def volume(x): # NOTE to be minimized
a = a_value*25.4/1000 # [m] along x
b = b_value*25.4/1000 # [m] along y
stack, thicknesses = layup_thick_from_x(x)
return sum(thicknesses)*a*b # [m^3]
Defining constraint functions#
[82]:
def calc_constr_buckling(x): # NOTE feasible when <= 0
stack, thicknesses = layup_thick_from_x(x)
prop = laminated_plate(stack=stack, plyts=thicknesses, laminaprop=laminaprop)
lambda_cb = calc_buckling_analytical(prop)
return 1 - lambda_cb/target_lambda
def calc_constr_failure(x): # NOTE feasible when <= 0
stack, thicknesses = layup_thick_from_x(x)
prop = laminated_plate(stack=stack, plyts=thicknesses, laminaprop=laminaprop)
lambda_cs = calc_failure_load_Haftka(prop)
return 1 - lambda_cs/target_lambda
test = [2, 1, 1, 1, 1, 0, 0, 1, 0, 0, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1]
stack, thicknesses = layup_thick_from_x(test)
print(stack.tolist())
print(stack_ref_half)
print(len(stack))
prop = laminated_plate(stack=stack, plyts=thicknesses, laminaprop=laminaprop)
lambda_cb = calc_buckling_analytical(prop)
print(lambda_cb)
[90, 90, 45, -45, 45, -45, 45, -45, 45, -45, 0, 0, 0, 0, 45, -45, 0, 0, 0, 0, 45, -45, 0, 0, 0, 0, -45, 45, 0, 0, 0, 0, -45, 45, 0, 0, 0, 0, -45, 45, -45, 45, -45, 45, -45, 45, 90, 90]
[90, 90, 45, -45, 45, -45, 45, -45, 45, -45, 0, 0, 0, 0, 45, -45, 0, 0, 0, 0, 45, -45, 0, 0]
48
14167.523183899806
Global optimization using pymoo#
Based on the tutorials:
https://pymoo.org/problems/definition.html
https://pymoo.org/customization/discrete.html?highlight=integerrandomsampling
https://www.pymoo.org/constraints/index.html
Here we define the optimization problem:
[83]:
from pymoo.core.problem import ElementwiseProblem
from pymoo.optimize import minimize
from pymoo.algorithms.soo.nonconvex.ga import GA
from pymoo.operators.sampling.rnd import IntegerRandomSampling
from pymoo.operators.crossover.sbx import SBX
from pymoo.operators.mutation.pm import PM
from pymoo.operators.repair.rounding import RoundingRepair
n_plies = 48
n_var = 2*(n_plies//4)
# define permutation problem
class Problem(ElementwiseProblem):
def __init__(self):
super().__init__(n_var = n_var,
n_obj = 1,
xl = [0]*n_var, # NOTE variables always changing from 0 to xu
xu = [len(available_angles)-1]*(n_var//2) + [1]*(n_var//2),
n_ieq_constr = 2,
vtype = int)
def _evaluate(self, x, out, *args, **kwargs):
x = np.asarray(x, dtype=int)
layup, plyts = layup_thick_from_x(x)
out['F'] = volume(x) # objective function
out['G'] = [calc_constr_buckling(x), calc_constr_failure(x)]
problem = Problem()
Defining the optimization method, here using a genetic algorithm:
[87]:
method = algorithm = GA(
pop_size=20,
sampling=IntegerRandomSampling(),
crossover=SBX(prob=1.0, eta=3.0, vtype=float, repair=RoundingRepair()),
mutation=PM(prob=1.0, eta=3.0, vtype=float, repair=RoundingRepair()),
eliminate_duplicates=True)
Running the optimization:
[88]:
res = minimize(problem,
method,
termination=('n_gen', 40),
seed=1,
save_history=True
)
print("Best solution found: %s" % res.X)
print("Function value: %s" % res.F)
print("Constraint violation: %s" % res.CV)
/var/folders/x7/jww7lc053bv10s9m981_324dmy7lhy/T/ipykernel_57190/2488278243.py:25: RuntimeWarning: divide by zero encountered in divide
gamma_12_allowable/abs(gamma_i_12) - 1
Best solution found: [1 1 2 2 0 0 2 1 0 0 0 0 1 1 1 0 1 1 1 1 1 1 1 1]
Function value: [0.00036052]
Constraint violation: [0.]
Checking the optimization result:
[90]:
stack_opt, thicknesses_opt = layup_thick_from_x(res.X)
prop_opt = laminated_plate(stack=stack_opt, plyts=thicknesses_opt, laminaprop=laminaprop)
ms_buckling_opt = calc_buckling_analytical(prop_opt)/target_lambda - 1
ms_failure_opt = calc_failure_load_Haftka(prop_opt)/target_lambda - 1
print('num of layers', len(stack_opt))
print('volume', volume(res.X))
print('ms_buckling_opt', ms_buckling_opt)
print('ms_failure_opt', ms_failure_opt)
print('optimal layup', stack_opt)
num of layers 44
volume 0.0003605154080000002
ms_buckling_opt 0.023418794562162004
ms_failure_opt 0.4356104684218065
optimal layup [ 45 -45 45 -45 90 90 0 0 0 0 90 90 45 -45 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 -45 45 90 90 0 0
0 0 90 90 -45 45 -45 45]
/var/folders/x7/jww7lc053bv10s9m981_324dmy7lhy/T/ipykernel_57190/2488278243.py:25: RuntimeWarning: divide by zero encountered in divide
gamma_12_allowable/abs(gamma_i_12) - 1
[ ]: