Revisiting the Beam Problem - Minimizing Stress with KS Constraints and BSplines#
The following example shows the optimization of a simple beam with rectangular cross section. The beam has been subdivided into elements along the length of the beam, and one end is fixed. The goal is to minimize the volume (and hence the mass of the homogeneous beam) by varying the thickness in each element without exceeding a maximum stress constraint while the beam is subject to multiple load cases, each one being a distributed force load that varies sinusoidally along the span.
Constraining the bending stress on each element leads to a more computationally expensive derivative calculation, so we
will use the KSFunction to reduce the stress vector for each load case to a single scalar value. To do so, we also need
to insert an ExecComp component that converts the stress into a form where a negative value means it is satisfied, and
a positive value means it is violated.
The problem presented here is also an example of a multi-point implementation, where we create a separate instance of the parts of the calculation that are impacted by different load cases. This enables our model to take advantage of multiple processors when run in parallel.
If we allow the optimizer to vary the thickness of each element, then we have a design variable vector that is as wide as the
number of elements in the model. This may perform poorly if we have a large number of elements in the beam. If we assume that
the optimal beam thickness is going to have a smooth continuous variation over the length of the beam, then it is a good
candidate for using an interpolation component like SplineComp to reduce the number of design variables we need.
For this example, we have 25 elements, but can reduce that to 5 control points for the optimizer’s design variables by including the SplineComp.
The code for the top system is this:
LocalStiffnessMatrixComp class definition
class LocalStiffnessMatrixComp(om.ExplicitComponent):
def initialize(self):
self.options.declare('num_elements', types=int)
self.options.declare('E')
self.options.declare('L')
def setup(self):
num_elements = self.options['num_elements']
E = self.options['E']
L = self.options['L']
self.add_input('I', shape=num_elements)
self.add_output('K_local', shape=(num_elements, 4, 4))
L0 = L / num_elements
coeffs = np.empty((4, 4))
coeffs[0, :] = [12, 6 * L0, -12, 6 * L0]
coeffs[1, :] = [6 * L0, 4 * L0 ** 2, -6 * L0, 2 * L0 ** 2]
coeffs[2, :] = [-12, -6 * L0, 12, -6 * L0]
coeffs[3, :] = [6 * L0, 2 * L0 ** 2, -6 * L0, 4 * L0 ** 2]
coeffs *= E / L0 ** 3
self.mtx = np.zeros((num_elements, 4, 4, num_elements))
for ind in range(num_elements):
self.mtx[ind, :, :, ind] = coeffs
self.declare_partials('K_local', 'I',
val=self.mtx.reshape(16 * num_elements, num_elements))
def compute(self, inputs, outputs):
outputs['K_local'] = 0
for ind in range(self.options['num_elements']):
outputs['K_local'][ind, :, :] = self.mtx[ind, :, :, ind] * inputs['I'][ind]
MultiStatesComp class definition
class MultiStatesComp(om.ImplicitComponent):
def initialize(self):
self.options.declare('num_elements', types=int)
self.options.declare('force_vector', types=np.ndarray)
self.options.declare('num_rhs', types=int)
def setup(self):
num_elements = self.options['num_elements']
num_nodes = num_elements + 1
size = 2 * num_nodes + 2
num_rhs = self.options['num_rhs']
self.add_input('K_local', shape=(num_elements, 4, 4))
for j in range(num_rhs):
self.add_output('d_%d' % j, shape=size)
cols = np.arange(16*num_elements)
rows = np.repeat(np.arange(4), 4)
rows = np.tile(rows, num_elements) + np.repeat(np.arange(num_elements), 16) * 2
for j in range(num_rhs):
disp = 'd_%d' % j
self.declare_partials(disp, 'K_local', rows=rows, cols=cols)
self.declare_partials(disp, disp)
def apply_nonlinear(self, inputs, outputs, residuals):
num_rhs = self.options['num_rhs']
self.K = self.assemble_CSC_K(inputs)
for j in range(num_rhs):
force_vector = np.concatenate([self.options['force_vector'][:, j], np.zeros(2)])
residuals['d_%d' % j] = self.K.dot(outputs['d_%d' % j]) - force_vector
def solve_nonlinear(self, inputs, outputs):
num_rhs = self.options['num_rhs']
self.K = self.assemble_CSC_K(inputs)
self.lu = splu(self.K)
for j in range(num_rhs):
force_vector = np.concatenate([self.options['force_vector'][:, j], np.zeros(2)])
outputs['d_%d' % j] = self.lu.solve(force_vector)
def linearize(self, inputs, outputs, jacobian):
num_rhs = self.options['num_rhs']
num_elements = self.options['num_elements']
self.K = self.assemble_CSC_K(inputs)
self.lu = splu(self.K)
i_elem = np.tile(np.arange(4), 4)
i_d = np.tile(i_elem, num_elements) + np.repeat(np.arange(num_elements), 16) * 2
K_dense = self.K.toarray()
for j in range(num_rhs):
disp = 'd_%d' % j
jacobian[disp, 'K_local'] = outputs[disp][i_d]
jacobian[disp, disp] = K_dense
def solve_linear(self, d_outputs, d_residuals, mode):
num_rhs = self.options['num_rhs']
for j in range(num_rhs):
disp = 'd_%d' % j
if mode == 'fwd':
d_outputs[disp] = self.lu.solve(d_residuals[disp])
else:
d_residuals[disp] = self.lu.solve(d_outputs[disp])
def assemble_CSC_K(self, inputs):
"""
Assemble the stiffness matrix in sparse CSC format.
Returns
-------
ndarray
Stiffness matrix as dense ndarray.
"""
num_elements = self.options['num_elements']
num_nodes = num_elements + 1
num_entry = num_elements * 12 + 4
ndim = num_entry + 4
data = np.zeros((ndim, ), dtype=inputs._get_data().dtype)
cols = np.empty((ndim, ))
rows = np.empty((ndim, ))
# First element.
data[:16] = inputs['K_local'][0, :, :].flat
cols[:16] = np.tile(np.arange(4), 4)
rows[:16] = np.repeat(np.arange(4), 4)
j = 16
for ind in range(1, num_elements):
ind1 = 2 * ind
K = inputs['K_local'][ind, :, :]
# NW quadrant gets summed with previous connected element.
data[j-6:j-4] += K[0, :2]
data[j-2:j] += K[1, :2]
# NE quadrant
data[j:j+4] = K[:2, 2:].flat
rows[j:j+4] = np.array([ind1, ind1, ind1 + 1, ind1 + 1])
cols[j:j+4] = np.array([ind1 + 2, ind1 + 3, ind1 + 2, ind1 + 3])
# SE and SW quadrants together
data[j+4:j+12] = K[2:, :].flat
rows[j+4:j+12] = np.repeat(np.arange(ind1 + 2, ind1 + 4), 4)
cols[j+4:j+12] = np.tile(np.arange(ind1, ind1 + 4), 2)
j += 12
data[-4:] = 1.0
rows[-4] = 2 * num_nodes
rows[-3] = 2 * num_nodes + 1
rows[-2] = 0.0
rows[-1] = 1.0
cols[-4] = 0.0
cols[-3] = 1.0
cols[-2] = 2 * num_nodes
cols[-1] = 2 * num_nodes + 1
n_K = 2 * num_nodes + 2
return coo_matrix((data, (rows, cols)), shape=(n_K, n_K)).tocsc()
MomentOfInertiaComp class definition
class MomentOfInertiaComp(om.ExplicitComponent):
def initialize(self):
self.options.declare('num_elements', types=int)
self.options.declare('b')
def setup(self):
num_elements = self.options['num_elements']
self.add_input('h', shape=num_elements)
self.add_output('I', shape=num_elements)
def setup_partials(self):
rows = cols = np.arange(self.options['num_elements'])
self.declare_partials('I', 'h', rows=rows, cols=cols)
def compute(self, inputs, outputs):
outputs['I'] = 1./12. * self.options['b'] * inputs['h'] ** 3
def compute_partials(self, inputs, partials):
partials['I', 'h'] = 1./4. * self.options['b'] * inputs['h'] ** 2
MultiStressComp class definition
class MultiStressComp(om.ExplicitComponent):
def initialize(self):
self.options.declare('num_elements', types=int)
self.options.declare('num_rhs', types=int)
self.options.declare('E')
def setup(self):
num_elements = self.options['num_elements']
num_nodes = num_elements + 1
num_rhs = self.options['num_rhs']
self.add_input('h', shape=num_elements)
for j in range(num_rhs):
self.add_input('displacements_%d' % j, shape=2 * num_nodes)
self.add_output('stress_%d' % j, shape=num_elements)
self.declare_partials(of='stress_%d' % j, wrt='displacements_%d' % j)
self.declare_partials(of='stress_%d' % j, wrt='h')
def compute(self, inputs, outputs):
num_rhs = self.options['num_rhs']
tk = inputs['h'] * 0.5
E = self.options['E']
for j in range(num_rhs):
ang = inputs['displacements_%d' % j][1::2]
d_ang = ang[1:] - ang[0:-1]
outputs['stress_%d' % j] = tk * E * d_ang
def compute_partials(self, inputs, partials):
num_elements = self.options['num_elements']
num_nodes = num_elements + 1
num_states = num_nodes * 2
num_rhs = self.options['num_rhs']
tk = inputs['h'] * 0.5
E = self.options['E']
for j in range(num_rhs):
ang = inputs['displacements_%d' % j][1::2]
d_ang = ang[1:] - ang[0:-1]
partials['stress_%d' % j, 'h'] = np.diag(0.5 * E * d_ang)
J = np.zeros((num_elements, num_states))
J[:, 1:-1:2] = -np.diag(tk * E)
J[:, 3::2] += np.diag(tk * E)
partials['stress_%d' % j, 'displacements_%d' % j] = J
VolumeComp class definition
class VolumeComp(om.ExplicitComponent):
def initialize(self):
self.options.declare('num_elements', types=int)
self.options.declare('b', default=1.)
self.options.declare('L')
def setup(self):
num_elements = self.options['num_elements']
b = self.options['b']
L = self.options['L']
L0 = L / num_elements
self.add_input('h', shape=num_elements)
self.add_output('volume')
self.declare_partials('volume', 'h', val=b * L0)
def compute(self, inputs, outputs):
L0 = self.options['L'] / self.options['num_elements']
outputs['volume'] = np.sum(inputs['h'] * self.options['b'] * L0)
"""
This is a multipoint implementation of the beam optimization problem.
This version minimizes volume while satisfying a max bending stress constraint in each element
for each loadcase.
"""
import numpy as np
import openmdao.api as om
from openmdao.test_suite.test_examples.beam_optimization.components.local_stiffness_matrix_comp import LocalStiffnessMatrixComp
from openmdao.test_suite.test_examples.beam_optimization.components.moment_comp import MomentOfInertiaComp
from openmdao.test_suite.test_examples.beam_optimization.components.multi_stress_comp import MultiStressComp
from openmdao.test_suite.test_examples.beam_optimization.components.volume_comp import VolumeComp
from openmdao.utils.spline_distributions import sine_distribution
def divide_cases(ncases, nprocs):
"""
Divide up load cases among available procs.
Parameters
----------
ncases : int
Number of load cases.
nprocs : int
Number of processors.
Returns
-------
list of list of int
Integer case numbers for each proc.
"""
data = []
for j in range(nprocs):
data.append([])
wrap = 0
for j in range(ncases):
idx = j - wrap
if idx >= nprocs:
idx = 0
wrap = j
data[idx].append(j)
return data
class MultipointBeamGroup(om.Group):
"""
System setup for minimization of volume (i.e., mass) subject to KS aggregated bending stress constraints.
"""
def initialize(self):
self.options.declare('E')
self.options.declare('L')
self.options.declare('b')
self.options.declare('volume')
self.options.declare('max_bending')
self.options.declare('num_elements', 5)
self.options.declare('num_cp', 50)
self.options.declare('num_load_cases', 1)
self.options.declare('parallel_derivs', False, types=bool, allow_none=True)
self.options.declare('ks_add_constraint', default=False, types=bool)
def setup(self):
E = self.options['E']
L = self.options['L']
b = self.options['b']
# volume = self.options['volume']
max_bending = self.options['max_bending']
num_elements = self.options['num_elements']
num_nodes = num_elements + 1
num_cp = self.options['num_cp']
num_load_cases = self.options['num_load_cases']
parallel_derivs = self.options['parallel_derivs']
x_interp = sine_distribution(num_elements)
comp = om.SplineComp(method='bsplines', num_cp=num_cp, x_interp_val=x_interp)
comp.add_spline(y_cp_name='h_cp', y_interp_name='h')
self.add_subsystem('interp', comp)
I_comp = MomentOfInertiaComp(num_elements=num_elements, b=b)
self.add_subsystem('I_comp', I_comp)
comp = LocalStiffnessMatrixComp(num_elements=num_elements, E=E, L=L)
self.add_subsystem('local_stiffness_matrix_comp', comp)
# Parallel Subsystem for load cases.
par = self.add_subsystem('parallel', om.ParallelGroup())
# Determine how to split cases up over the available procs.
nprocs = self.comm.size
divide = divide_cases(num_load_cases, nprocs)
for j, this_proc in enumerate(divide):
num_rhs = len(this_proc)
name = 'sub_%d' % j
sub = par.add_subsystem(name, om.Group())
# Load is a sinusoidal distributed force of varying spatial frequency.
force_vector = np.zeros((2 * num_nodes, num_rhs))
for i, k in enumerate(this_proc):
end = 1.5 * np.pi
if num_load_cases > 1:
end += k * 0.5 * np.pi / (num_load_cases - 1)
x = np.linspace(0, end, num_nodes)
f = - np.sin(x)
force_vector[0:-1:2, i] = f
comp = MultiStatesComp(num_elements=num_elements, force_vector=force_vector,
num_rhs=num_rhs)
sub.add_subsystem('states_comp', comp)
comp = MultiStressComp(num_elements=num_elements, E=E, num_rhs=num_rhs)
sub.add_subsystem('stress_comp', comp)
self.connect('local_stiffness_matrix_comp.K_local',
'parallel.%s.states_comp.K_local' % name)
for k in range(num_rhs):
sub.connect('states_comp.d_%d' % k,
'stress_comp.displacements_%d' % k,
src_indices=np.arange(2 *num_nodes))
if parallel_derivs:
color = 'red_%d' % k
else:
color = None
comp = om.KSComp(width=num_elements, upper=max_bending,
add_constraint=self.options['ks_add_constraint'],
parallel_deriv_color=color)
sub.add_subsystem('KS_%d' % k, comp)
sub.connect('stress_comp.stress_%d' % k,
'KS_%d.g' % k)
if not self.options['ks_add_constraint']:
sub.add_constraint('KS_%d.KS' % k, upper=0.0,
parallel_deriv_color=color)
comp = VolumeComp(num_elements=num_elements, b=b, L=L)
self.add_subsystem('volume_comp', comp)
self.connect('interp.h', 'I_comp.h')
self.connect('interp.h', 'volume_comp.h')
self.connect('I_comp.I', 'local_stiffness_matrix_comp.I')
self.add_design_var('interp.h_cp', lower=1e-2, upper=10.)
self.add_objective('volume_comp.volume')
Next we run the model, and choose ScipyOptimizeDriver SLSQP to be our optimizer. At the conclusion of optimization, we print out the design variable, which is the thickness for each element.
from openmdao.test_suite.test_examples.beam_optimization.components.multi_states_comp import MultiStatesComp
E = 1.
L = 1.
b = 0.1
volume = 0.01
max_bending = 100.0
num_cp = 5
num_elements = 25
num_load_cases = 2
model = MultipointBeamGroup(E=E, L=L, b=b, volume=volume, max_bending = max_bending,
num_elements=num_elements, num_cp=num_cp,
num_load_cases=num_load_cases)
prob = om.Problem(model=model)
prob.driver = om.ScipyOptimizeDriver()
prob.driver.options['optimizer'] = 'SLSQP'
prob.driver.options['tol'] = 1e-9
prob.driver.options['disp'] = True
prob.setup(mode='rev')
prob.run_driver()
print(prob['interp.h'][0])
Optimization terminated successfully (Exit mode 0)
Current function value: 0.03494890159976332
Iterations: 40
Function evaluations: 117
Gradient evaluations: 40
Optimization Complete
-----------------------------------
[0.45632323 0.45612552 0.45543324 0.45397058 0.45134629 0.44714397
0.4410258 0.43283139 0.42265378 0.41087801 0.39817311 0.38543581
0.37369202 0.36342186 0.35289066 0.34008777 0.32362887 0.30300358
0.27867837 0.25204063 0.22519409 0.20063906 0.18088818 0.16807856
0.16364105]