{
"cells": [
{
"cell_type": "code",
"execution_count": 1,
"id": "509d0124",
"metadata": {
"execution": {
"iopub.execute_input": "2026-10-02T14:40:52.935820Z",
"iopub.status.busy": "2026-10-02T14:40:52.935509Z",
"iopub.status.idle": "2026-10-02T14:40:52.939370Z",
"shell.execute_reply": "2026-10-02T14:40:52.938917Z"
},
"hide_input": true,
"papermill": {
"duration": 0.00642,
"end_time": "2026-10-02T14:40:52.939796+00:00",
"exception": false,
"start_time": "2026-10-02T14:40:52.933376+00:00",
"status": "completed"
},
"tags": [
"active-ipynb",
"remove-input",
"remove-output"
]
},
"outputs": [],
"source": [
"try:\n",
" from openmdao.utils.notebook_utils import notebook_mode # noqa: F401\n",
"except ImportError:\n",
" !python -m pip install openmdao[notebooks]"
]
},
{
"cell_type": "markdown",
"id": "3ed15d15",
"metadata": {
"papermill": {
"duration": 0.001032,
"end_time": "2026-10-02T14:40:52.942229+00:00",
"exception": false,
"start_time": "2026-10-02T14:40:52.941197+00:00",
"status": "completed"
},
"tags": []
},
"source": [
"# Revisiting the Beam Problem - Minimizing Stress with KS Constraints and BSplines\n",
"\n",
"The following example shows the optimization of a simple beam with rectangular cross section. The beam has been subdivided into\n",
"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\n",
"homogeneous beam) by varying the thickness in each element without exceeding a maximum stress constraint while the beam is\n",
"subject to multiple load cases, each one being a distributed force load that varies sinusoidally along the span.\n",
"\n",
"Constraining the bending stress on each element leads to a more computationally expensive derivative calculation, so we\n",
"will use the `KSFunction` to reduce the stress vector for each load case to a single scalar value. To do so, we also need\n",
"to insert an `ExecComp` component that converts the stress into a form where a negative value means it is satisfied, and\n",
"a positive value means it is violated.\n",
"\n",
"The problem presented here is also an example of a multi-point implementation, where we create a separate instance of the\n",
"parts of the calculation that are impacted by different load cases. This enables our model to take advantage of multiple\n",
"processors when run in parallel.\n",
"\n",
"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\n",
"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\n",
"the optimal beam thickness is going to have a smooth continuous variation over the length of the beam, then it is a good\n",
"candidate for using an interpolation component like `SplineComp` to reduce the number of design variables we need.\n",
"\n",
"For this example, we have 25 elements, but can reduce that to 5 control points for the optimizer's design variables by\n",
"including the SplineComp.\n",
"\n",
"The code for the top system is this:"
]
},
{
"cell_type": "code",
"execution_count": 2,
"id": "1be9f90d",
"metadata": {
"execution": {
"iopub.execute_input": "2026-10-02T14:40:52.945082Z",
"iopub.status.busy": "2026-10-02T14:40:52.944835Z",
"iopub.status.idle": "2026-10-02T14:40:55.594690Z",
"shell.execute_reply": "2026-10-02T14:40:55.594026Z"
},
"papermill": {
"duration": 2.652057,
"end_time": "2026-10-02T14:40:55.595343+00:00",
"exception": false,
"start_time": "2026-10-02T14:40:52.943286+00:00",
"status": "completed"
},
"tags": [
"remove-input",
"remove-output"
]
},
"outputs": [
{
"data": {
"application/papermill.record/text/html": "class LocalStiffnessMatrixComp(om.ExplicitComponent):\n\n def initialize(self):\n self.options.declare('num_elements', types=int)\n self.options.declare('E')\n self.options.declare('L')\n\n def setup(self):\n num_elements = self.options['num_elements']\n E = self.options['E']\n L = self.options['L']\n\n self.add_input('I', shape=num_elements)\n self.add_output('K_local', shape=(num_elements, 4, 4))\n\n L0 = L / num_elements\n coeffs = np.empty((4, 4))\n coeffs[0, :] = [12, 6 * L0, -12, 6 * L0]\n coeffs[1, :] = [6 * L0, 4 * L0 ** 2, -6 * L0, 2 * L0 ** 2]\n coeffs[2, :] = [-12, -6 * L0, 12, -6 * L0]\n coeffs[3, :] = [6 * L0, 2 * L0 ** 2, -6 * L0, 4 * L0 ** 2]\n coeffs *= E / L0 ** 3\n\n self.mtx = np.zeros((num_elements, 4, 4, num_elements))\n for ind in range(num_elements):\n self.mtx[ind, :, :, ind] = coeffs\n\n self.declare_partials('K_local', 'I',\n val=self.mtx.reshape(16 * num_elements, num_elements))\n\n def compute(self, inputs, outputs):\n outputs['K_local'] = 0\n for ind in range(self.options['num_elements']):\n outputs['K_local'][ind, :, :] = self.mtx[ind, :, :, ind] * inputs['I'][ind]\nclass MultiStatesComp(om.ImplicitComponent):\n\n def initialize(self):\n self.options.declare('num_elements', types=int)\n self.options.declare('force_vector', types=np.ndarray)\n self.options.declare('num_rhs', types=int)\n\n def setup(self):\n num_elements = self.options['num_elements']\n num_nodes = num_elements + 1\n size = 2 * num_nodes + 2\n num_rhs = self.options['num_rhs']\n\n self.add_input('K_local', shape=(num_elements, 4, 4))\n for j in range(num_rhs):\n self.add_output('d_%d' % j, shape=size)\n\n cols = np.arange(16*num_elements)\n rows = np.repeat(np.arange(4), 4)\n rows = np.tile(rows, num_elements) + np.repeat(np.arange(num_elements), 16) * 2\n\n for j in range(num_rhs):\n disp = 'd_%d' % j\n\n self.declare_partials(disp, 'K_local', rows=rows, cols=cols)\n self.declare_partials(disp, disp)\n\n def apply_nonlinear(self, inputs, outputs, residuals):\n num_rhs = self.options['num_rhs']\n\n self.K = self.assemble_CSC_K(inputs)\n for j in range(num_rhs):\n force_vector = np.concatenate([self.options['force_vector'][:, j], np.zeros(2)])\n residuals['d_%d' % j] = self.K.dot(outputs['d_%d' % j]) - force_vector\n\n def solve_nonlinear(self, inputs, outputs):\n num_rhs = self.options['num_rhs']\n\n self.K = self.assemble_CSC_K(inputs)\n self.lu = splu(self.K)\n\n for j in range(num_rhs):\n\n force_vector = np.concatenate([self.options['force_vector'][:, j], np.zeros(2)])\n outputs['d_%d' % j] = self.lu.solve(force_vector)\n\n def linearize(self, inputs, outputs, jacobian):\n num_rhs = self.options['num_rhs']\n num_elements = self.options['num_elements']\n\n self.K = self.assemble_CSC_K(inputs)\n self.lu = splu(self.K)\n\n i_elem = np.tile(np.arange(4), 4)\n i_d = np.tile(i_elem, num_elements) + np.repeat(np.arange(num_elements), 16) * 2\n\n K_dense = self.K.toarray()\n\n for j in range(num_rhs):\n disp = 'd_%d' % j\n jacobian[disp, 'K_local'] = outputs[disp][i_d]\n jacobian[disp, disp] = K_dense\n\n def solve_linear(self, d_outputs, d_residuals, mode):\n num_rhs = self.options['num_rhs']\n\n for j in range(num_rhs):\n disp = 'd_%d' % j\n\n if mode == 'fwd':\n d_outputs[disp] = self.lu.solve(d_residuals[disp])\n else:\n d_residuals[disp] = self.lu.solve(d_outputs[disp])\n\n def assemble_CSC_K(self, inputs):\n """\n Assemble the stiffness matrix in sparse CSC format.\n\n Returns\n -------\n ndarray\n Stiffness matrix as dense ndarray.\n """\n num_elements = self.options['num_elements']\n num_nodes = num_elements + 1\n num_entry = num_elements * 12 + 4\n ndim = num_entry + 4\n\n data = np.zeros((ndim, ), dtype=inputs._get_data().dtype)\n cols = np.empty((ndim, ))\n rows = np.empty((ndim, ))\n\n # First element.\n data[:16] = inputs['K_local'][0, :, :].flat\n cols[:16] = np.tile(np.arange(4), 4)\n rows[:16] = np.repeat(np.arange(4), 4)\n\n j = 16\n for ind in range(1, num_elements):\n ind1 = 2 * ind\n K = inputs['K_local'][ind, :, :]\n\n # NW quadrant gets summed with previous connected element.\n data[j-6:j-4] += K[0, :2]\n data[j-2:j] += K[1, :2]\n\n # NE quadrant\n data[j:j+4] = K[:2, 2:].flat\n rows[j:j+4] = np.array([ind1, ind1, ind1 + 1, ind1 + 1])\n cols[j:j+4] = np.array([ind1 + 2, ind1 + 3, ind1 + 2, ind1 + 3])\n\n # SE and SW quadrants together\n data[j+4:j+12] = K[2:, :].flat\n rows[j+4:j+12] = np.repeat(np.arange(ind1 + 2, ind1 + 4), 4)\n cols[j+4:j+12] = np.tile(np.arange(ind1, ind1 + 4), 2)\n\n j += 12\n\n data[-4:] = 1.0\n rows[-4] = 2 * num_nodes\n rows[-3] = 2 * num_nodes + 1\n rows[-2] = 0.0\n rows[-1] = 1.0\n cols[-4] = 0.0\n cols[-3] = 1.0\n cols[-2] = 2 * num_nodes\n cols[-1] = 2 * num_nodes + 1\n\n n_K = 2 * num_nodes + 2\n return coo_matrix((data, (rows, cols)), shape=(n_K, n_K)).tocsc()\nclass MomentOfInertiaComp(om.ExplicitComponent):\n\n def initialize(self):\n self.options.declare('num_elements', types=int)\n self.options.declare('b')\n\n def setup(self):\n num_elements = self.options['num_elements']\n\n self.add_input('h', shape=num_elements)\n self.add_output('I', shape=num_elements)\n\n def setup_partials(self):\n rows = cols = np.arange(self.options['num_elements'])\n self.declare_partials('I', 'h', rows=rows, cols=cols)\n\n def compute(self, inputs, outputs):\n outputs['I'] = 1./12. * self.options['b'] * inputs['h'] ** 3\n\n def compute_partials(self, inputs, partials):\n partials['I', 'h'] = 1./4. * self.options['b'] * inputs['h'] ** 2\nclass MultiStressComp(om.ExplicitComponent):\n\n def initialize(self):\n self.options.declare('num_elements', types=int)\n self.options.declare('num_rhs', types=int)\n self.options.declare('E')\n\n def setup(self):\n num_elements = self.options['num_elements']\n num_nodes = num_elements + 1\n num_rhs = self.options['num_rhs']\n\n self.add_input('h', shape=num_elements)\n\n for j in range(num_rhs):\n self.add_input('displacements_%d' % j, shape=2 * num_nodes)\n self.add_output('stress_%d' % j, shape=num_elements)\n\n self.declare_partials(of='stress_%d' % j, wrt='displacements_%d' % j)\n self.declare_partials(of='stress_%d' % j, wrt='h')\n\n def compute(self, inputs, outputs):\n num_rhs = self.options['num_rhs']\n tk = inputs['h'] * 0.5\n E = self.options['E']\n\n for j in range(num_rhs):\n ang = inputs['displacements_%d' % j][1::2]\n d_ang = ang[1:] - ang[0:-1]\n outputs['stress_%d' % j] = tk * E * d_ang\n\n def compute_partials(self, inputs, partials):\n num_elements = self.options['num_elements']\n num_nodes = num_elements + 1\n num_states = num_nodes * 2\n num_rhs = self.options['num_rhs']\n tk = inputs['h'] * 0.5\n E = self.options['E']\n\n for j in range(num_rhs):\n ang = inputs['displacements_%d' % j][1::2]\n d_ang = ang[1:] - ang[0:-1]\n partials['stress_%d' % j, 'h'] = np.diag(0.5 * E * d_ang)\n\n J = np.zeros((num_elements, num_states))\n J[:, 1:-1:2] = -np.diag(tk * E)\n J[:, 3::2] += np.diag(tk * E)\n partials['stress_%d' % j, 'displacements_%d' % j] = J\nclass VolumeComp(om.ExplicitComponent):\n\n def initialize(self):\n self.options.declare('num_elements', types=int)\n self.options.declare('b', default=1.)\n self.options.declare('L')\n\n def setup(self):\n num_elements = self.options['num_elements']\n b = self.options['b']\n L = self.options['L']\n L0 = L / num_elements\n\n self.add_input('h', shape=num_elements)\n self.add_output('volume')\n\n self.declare_partials('volume', 'h', val=b * L0)\n\n def compute(self, inputs, outputs):\n L0 = self.options['L'] / self.options['num_elements']\n\n outputs['volume'] = np.sum(inputs['h'] * self.options['b'] * L0)\n