Tutorial: Lightweight optimisation using a gradient-based lamination parameter approach#
Date: 4 February 2026
Author: Saullo G. P. Castro
Cite this tutorial as:
Castro, SGP. Methods for analysis and design of composites (Version 0.8.4) [Computer software]. 2024. https://doi.org/10.5281/zenodo.2871782
Installing required modules#
[13]:
!python -m pip install numpy scipy composites > tmp.txt
DEPRECATION: Loading egg at /Users/saullogiovanip/miniconda3/lib/python3.13/site-packages/panels-0.5.0-py3.13-macosx-11.1-arm64.egg is deprecated. pip 25.1 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.13/site-packages/tuduam-2026.3-py3.13.egg is deprecated. pip 25.1 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
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\):
[14]:
# 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:
[15]:
import numpy as np
# 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)
[16]:
def get_Q_matrix(theta_deg):
"""
Calculate the transformed reduced stiffness matrix Q_bar for a given angle.
"""
nu21 = nu12 * E22 / E11
denom = 1 - nu12 * nu21
Q11 = E11 / denom
Q22 = E22 / denom
Q12 = nu12 * E22 / denom
Q66 = G12
c = np.cos(np.radians(theta_deg))
s = np.sin(np.radians(theta_deg))
c2 = c**2; s2 = s**2; c4 = c**4; s4 = s**4
Q_bar_11 = Q11*c4 + 2*(Q12 + 2*Q66)*s2*c2 + Q22*s4
Q_bar_22 = Q11*s4 + 2*(Q12 + 2*Q66)*s2*c2 + Q22*c4
Q_bar_12 = (Q11 + Q22 - 4*Q66)*s2*c2 + Q12*(s4 + c4)
Q_bar_66 = (Q11 + Q22 - 2*Q12 - 2*Q66)*s2*c2 + Q66*(s4 + c4)
Q_bar_44 = G23 * c2 + G13 * s2
Q_bar_55 = G13 * c2 + G23 * s2
Q_bar_45 = (G13 - G23) * c * s
Q_bar = np.array([
[Q_bar_11, Q_bar_12, 0, 0, 0, 0],
[Q_bar_12, Q_bar_22, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0],
[0, 0, 0, Q_bar_44, Q_bar_45, 0],
[0, 0, 0, Q_bar_45, Q_bar_55, 0],
[0, 0, 0, 0, 0, Q_bar_66]
])
return Q_bar
[17]:
class LaminationParameterProp:
def __init__(self, h, xiD1, xiD3):
# Pre-calculate Base Q matrices
self.Q0 = get_Q_matrix(0)
Q11 = self.Q0[0,0]
Q22 = self.Q0[1,1]
Q12 = self.Q0[0,1]
Q66 = self.Q0[5,5]
self.xiD1 = xiD1
self.xiD2 = None
self.xiD3 = xiD3
self.xiD4 = None
self.U1 = (3*Q11 + 3*Q22 + 2*Q12 + 4*Q66)/8
self.U2 = (Q11 - Q22)/2
self.U3 = (Q11 + Q22 - 2*Q12 - 4*Q66)/8
self.U4 = (Q11 + Q22 + 6*Q12 - 4*Q66)/8
self.U5 = (Q11 + Q22 - 2*Q12 + 4*Q66)/8
self.h = h
self.D11 = h**3/12*(self.U1 + xiD1*self.U2 + xiD3*self.U3)
self.D22 = h**3/12*(self.U1 - xiD1*self.U2 + xiD3*self.U3)
self.D12 = h**3/12*(-xiD3*self.U3 + self.U4)
self.D66 = h**3/12*(-xiD3*self.U3 + self.U5)
Analytical equation to calculate buckling#
[18]:
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.D11
D12 = prop.D12
D22 = prop.D22
D66 = prop.D66
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
Analytical buckling gradients#
[19]:
def calc_buckling_analytical_jac(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.D11
D12 = prop.D12
D22 = prop.D22
D66 = prop.D66
lambda_b_min = 1e30
dlambda_b_dh = None
dlambda_b_dxiD1 = None
dlambda_b_dxiD3 = None
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)
)
lambdas = [lambda_b_min, lambda_b]
argmin = np.argmin(lambdas)
lambda_b_min = lambdas[argmin]
den = (m/a**2)*Nxx + (n/b)**2*Nyy
dlambda_b_dD11 = pi**2*(m/a)**4/den
dlambda_b_dD12 = 2*pi**2*(m/a)**2*(n/b)**2/den
dlambda_b_dD22 = pi**2*(n/b)**4/den
dlambda_b_dD66 = 4*pi**2*(m/a)**2*(n/b)**2/den
# Derivatives of D11
dD11_dh = 3 * prop.h**2 / 12 * (prop.U1 + prop.xiD1 * prop.U2 + prop.xiD3 * prop.U3)
dD11_dxiD1 = prop.h**3 / 12 * prop.U2
#dD11_dxiD2 = 0
dD11_dxiD3 = prop.h**3 / 12 * prop.U3
#dD11_dxiD4 = 0
# Derivatives of D22
dD22_dh = 3 * prop.h**2 / 12 * (prop.U1 - prop.xiD1 * prop.U2 + prop.xiD3 * prop.U3)
dD22_dxiD1 = -prop.h**3 / 12 * prop.U2
#dD22_dxiD2 = 0
dD22_dxiD3 = prop.h**3 / 12 * prop.U3
#dD22_dxiD4 = 0
# Derivatives of D12
dD12_dh = 0
dD12_dxiD1 = 0
#dD12_dxiD2 = 0
dD12_dxiD3 = -prop.h**3 / 12 * prop.U3
#dD12_dxiD4 = 0
# Derivatives of D66
dD66_dh = 0
dD66_dxiD1 = 0
dD66_dxiD2 = 0
dD66_dxiD3 = -prop.h**3 / 12 * prop.U3
dD66_dxiD4 = 0
dlambda_b_dh_new = (dlambda_b_dD11 * dD11_dh
+ dlambda_b_dD12 * dD12_dh
+ dlambda_b_dD22 * dD22_dh
+ dlambda_b_dD66 * dD66_dh)
dlambda_b_dxiD1_new = (dlambda_b_dD11 * dD11_dxiD1
+ dlambda_b_dD22 * dD22_dxiD1
+ dlambda_b_dD12 * dD12_dxiD1
+ dlambda_b_dD66 * dD66_dxiD1)
dlambda_b_dxiD3_new = (dlambda_b_dD11 * dD11_dxiD3
+ dlambda_b_dD22 * dD22_dxiD3
+ dlambda_b_dD12 * dD12_dxiD3
+ dlambda_b_dD66 * dD66_dxiD3)
dlambda_b_dh = [dlambda_b_dh, dlambda_b_dh_new][argmin]
dlambda_b_dxiD1 = [dlambda_b_dxiD1, dlambda_b_dxiD1_new][argmin]
dlambda_b_dxiD3 = [dlambda_b_dxiD3, dlambda_b_dxiD3_new][argmin]
dlambda_b_dx = [dlambda_b_dh, dlambda_b_dxiD1, dlambda_b_dxiD3]
return dlambda_b_dx
Verifying constraint functions#
Compare with column 2 of Table 3 below, first row.
[20]:
from composites import laminated_plate
h = ply_thickness*len(stack_ref)
lam = laminated_plate(plyt=ply_thickness, stack=stack_ref, laminaprop=laminaprop)
lp = lam.calc_lamination_parameters()
prop_ref = LaminationParameterProp(*[h, lp.xiD1, lp.xiD3])
lambda_cb_ref_analytical = calc_buckling_analytical(prop_ref)
print('lambda_cb_ref_analytical', lambda_cb_ref_analytical)
lambda_cb_ref_analytical 14167.52318389982
Defining design load#
[21]:
target_lambda = 10000
Defining objective function#
[22]:
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
h = x[0]
return h*a*b # [m^3]
Defining constraint functions#
[23]:
def calc_constr_buckling(x): # NOTE feasible when >= 0
prop = LaminationParameterProp(*x)
lambda_b = calc_buckling_analytical(prop)
return lambda_b/target_lambda - 1
def calc_constr_buckling_jac(x):
prop = LaminationParameterProp(*x)
dlambda_b_dx = calc_buckling_analytical_jac(prop)
return np.asarray(dlambda_b_dx)/target_lambda
Optimisation using SLSQP#
Here we define the optimization problem:
[24]:
from scipy.optimize import minimize
from scipy.optimize import Bounds
n_var = 3
h_min = ply_thickness
h_max = 48*ply_thickness
x0 = [(h_min + h_max) / 2.0, 0, 0]
bounds = Bounds([h_min, -1, -1], [h_max, +1, +1])
# NOTE inequality constraints: fun(x) >= 0
for use_gradient in [False, True]:
if use_gradient:
constraints = [{'type': 'ineq', 'fun': calc_constr_buckling, 'jac': calc_constr_buckling_jac},
]
else:
constraints = [{'type': 'ineq', 'fun': calc_constr_buckling},
]
print()
print("Starting SLSQP Optimization with Scipy, gradient=%s" % str(use_gradient))
print()
res = minimize(
volume,
x0,
method='SLSQP',
bounds=bounds,
constraints=constraints,
options={'ftol': 1e-9, 'disp': True, 'maxiter': 200},
jac="2-point"
)
print(" Best solution found: %s" % res.x)
print(" Checking best solution")
prop_opt = LaminationParameterProp(*res.x)
ms_buckling_opt = calc_buckling_analytical(prop_opt)/target_lambda - 1
print(' volume', volume(res.x))
print(' ms_buckling_opt', ms_buckling_opt)
lambda_cb_opt_analytical = calc_buckling_analytical(prop_opt)
print(' lambda_cb_opt_analytical', lambda_cb_opt_analytical)
Starting SLSQP Optimization with Scipy, gradient=False
Optimization terminated successfully (Exit mode 0)
Current function value: 0.00038613887522723374
Iterations: 5
Function evaluations: 23
Gradient evaluations: 5
Best solution found: [ 0.00598516 -0.82031951 -0.80568312]
Checking best solution
volume 0.00038613887522723374
ms_buckling_opt 1.1040242053894644e-09
lambda_cb_opt_analytical 10000.000011040242
Starting SLSQP Optimization with Scipy, gradient=True
Optimization terminated successfully (Exit mode 0)
Current function value: 0.00034799813181355646
Iterations: 15
Function evaluations: 112
Gradient evaluations: 15
Best solution found: [ 0.00539398 -0.64654288 -0.64782038]
Checking best solution
volume 0.00034799813181355646
ms_buckling_opt 4.441986687586663e-06
lambda_cb_opt_analytical 10000.044419866876
[ ]: