Mathematical Background ======================= .. sectionauthor:: Daniel Mohr .. currentmodule:: ttrr The software package :py:mod:`ttrr` emerged from the scientific research for [Mohr2011]_, [MohrSteinMatziesKnapek2012]_ and [MohrSteinMatziesKnapek2013]_. First, I started with MatLab http://www.mathworks.de/ in my PhD research, but the problems were too hard for this great software package. Furthermore, MatLab does not support the solving of mixed integer programs on its own. Therefore I was looking for something else: I found GLPK http://www.gnu.org/software/glpk/ and Python http://www.python.org/ --- both are great software! Because it was simple, I did the visualization with MatLab http://www.mathworks.de/ in my PhD-thesis [Mohr2011]_. To write :py:mod:`ttrr`, I do not want to use a commercial back end for the visualization or moreover point you to such a software. But if you like, you can do it on your own. .. contents:: General Concept of Topology Optimization of Truss ------------------------------------------------- .. index:: topology optimization of truss, ground structure Truss topology optimization has its root in [Michell1904]_. The analytically stated and known (cf. [Przemieniecki1968]_; [Marti2003]_; [BendsoeSigmund2004]_) linear problem of the topology optimization of a truss with respect to its volume (:index:`plastic design`) is given by the objective function to be minimized .. math:: :nowrap: \begin{align*} \min_{s,w}{l^T s} \end{align*} with the bar lengths :math:`l \in \R^n` subject to the equilibrium condition .. math:: :nowrap: \begin{align*} C w = f \end{align*} with the reduced :index:`geometry matrix` :math:`C \in \R^{m \times n}`, the inner bar forces :math:`w \in \R^n` and the reduced applied external loads :math:`f \in \R^m`. Further given are the constraints .. math:: :nowrap: \begin{align*} \sigma_- s \leq w \leq \sigma_+ s \end{align*} for the linear elasticity by --- for all bars equal --- :index:`stress limits` :math:`0 > \sigma_- \in \R` and :math:`0 < \sigma_+ \in \R` (e. g. yield points for pressure and tension of the material all bars are made of), and the box constraint .. math:: :nowrap: \begin{align*} 0 \leq s \leq s_{\max} \end{align*} for the cross sections of the bars :math:`s \in \R^n` with the given maximal bar cross section :math:`s_{\max} \in \R^n`. [#]_ Here :math:`w` is seen as a design variable; otherwise we have :math:`s(w)` as function and therefore no linear program. If more than one load (e. g. the set :math:`F` of loads) is to be considered, one ends up at the known multiple load case. But how is it performed? How can we put more than one load in the above linear program? It is impossible! But we can ask for a robust solution --- this means we consider the worst case: A solution for all :math:`f^{(i)} \in F` is searched. Now we really end up at the well known multiple load case: [#]_ .. math:: :nowrap: \begin{align*} \mbox{obj. func.: }&V(s) = l^T s \rightarrow \min_{s,w} \\ \mbox{s.t.: } & \forall f^{(i)} \in F: & \\ & \qquad C w^{(i)} = f^{(i)} \in \R^m \\ & \qquad \sigma_- s \leq w^{(i)} \leq \sigma_+ s \\ & 0 \leq s \leq s_{\max} \end{align*} All this is done by :py:func:`ttrr.ttrr`. .. [#] The maximal cross section is no limitation --- it's a feature! It provides the possibility to limit the cross section; but if you do not want this set a huge upper limit --- but an adequate limit helps the optimization algorithms. .. [#] Sorry for the compactness, but you can look at the literature all over the world. :-) Robust Topology Optimization of Truss ------------------------------------- .. index:: robust, robustness For :index:`robust optimization` see [Ben-TalElGhaouiNemirovski2009]_. The implementation in :py:func:`ttrr.ttrr_robust` follows [MohrSteinMatziesKnapek2012]_. This means if perturbed loads :math:`\tilde{F}` are given, they are composed by some perturbations of the unperturbed force :math:`f`. Due to the linear programming we can restrict ourselves to the extreme points of :math:`\tilde{F}`. If we start the other way round with an additional assumption, it becomes clear: Let :math:`f` be a given force with a perturbation of :math:`10 \mbox{\%}`. Then we are interested in a solution for all :math:`\tilde{f} \in \{\tilde{f} : 0.9 f \leq \tilde{f} \leq 1.1 f\}`. Obviously the extreme points are :math:`0.9 f` and :math:`1.1 f`. Therefore a solution for all loads in the set :math:`F := \{0.9 f ; 1.1 f\}` is searched. So robust optimization in topology optimization is the same as the multiple load case with appropriate load cases. For :math:`F` a polyhedron a finite number of extreme points is given. And therefore :py:func:`ttrr.ttrr` can be recycled. In addition :py:func:`ttrr.ttrr_robust` provides features to calculate the extreme points of simple sets of perturbed forces. Redundant Topology Optimization of Truss ---------------------------------------- The basic concept of topology optimization with regard to :index:`redundancy` is described in [Mohr2011]_ and in [MohrSteinMatziesKnapek2013]_. :index:`Redundant Topology Optimization of Truss` means we are looking for the best topology of a truss with a redundancy as to the falling out of some parts of the structure. The dream is a structure which is functional even if a significant damage occurs to an arbitrary set of bars. You should think of RAID [#]_. [PattersonGibsonKatz1987]_ For example, a :index:`RAID 1` is a bunch of discs and every disc does the work. This means every other one can fail. Following [Mohr2011]_ this is called :index:`"1-1/n"-redundancy` (c.f. [Shannon1948]_ [PattersonGibsonKatz1987]_). Transferred to structures this means we are looking for :math:`n` different structures in the same design space fulfilling the task to sustain the load(s). Let :math:`s_i \in \R^N` be the cross sections of the bars for the structure :math:`i`. The sum :math:`s = \sum_{i=1}^{n}{s_i}` describes the cross sections of the bars for the complete structure. The complete structure is composed of the :math:`n` structures :math:`s_i`. With the selecting variables :math:`x_i \in \{0;1\}^N` we can now formulate the problem as a mixed integer linear program: [#2]_ .. math:: :nowrap: \begin{align*} \mbox{obj. func.: }&V(s_1,s_2,\ldots,s_n) = \sum_{i=1}^{n}{l^T s_i} \rightarrow \min_{x,w,s} \\ \mbox{s.t.: } & \forall i=1,2,\ldots,n: \\ & \qquad s_i \leq s_{\max} x_i \in \R^N \\ & \qquad \forall k=1,2,\ldots,n_l: \\ & \qquad \qquad C w_i^{(k)} = f^{(k)} \in \R^m \\ & \qquad \qquad \sigma_- s_i \leq w_i^{(k)} \leq \sigma_+ s_i \\ & \sum_{i=1}^{n}{x_i} \leq 1 \\ & x \in \{0;1\}^{nN} ; \qquad 0 \leq s \end{align*} From :math:`x \in \{0;1\}^{nN}` follows :math:`0 \leq \sum_{i=1}^{n}{x_i}` and :math:`s \leq s_{\max}`. This is done by :py:func:`ttrr.ttrr_redundant` for *redundancy* = :math:`-n \leq -2`. A :index:`RAID 5` (c.f. [PattersonGibsonKatz1987]_) is more complex. It can be described as a simple code. To transfer this concept to structures we only need to know the behavior. A RAID 5 is a bunch of discs. Every single disc may fail and the others do the job together. Following [Mohr2011]_ this is called :index:`"1/n"-redundancy` (c.f. [Shannon1948]_ [PattersonGibsonKatz1987]_). Back to structures this means we are looking for :math:`n` different structures in the same design space. Every single structure can fall out and the other ones should sustain the load(s) together. Let :math:`s_i \in \R^N` be the cross sections of the bars for the structure :math:`i`. As before :math:`s = \sum_{i=1}^{n}{s_i}` describes the cross sections of the bars for the complete structure. The complete structure :math:`s` is composed of the :math:`n` structures :math:`s_i`. If a structure :math:`i` fall out, the resulting structure is :math:`\bar{s}_i = \sum_{j=1}^{n}{s_j} - s_i`. As mentioned before, this damaged structure should sustain the load(s). With the selecting variables :math:`x_i \in \{0;1\}^N` we can now formulate the problem as a mixed integer linear program: [#2]_ .. math:: :nowrap: \begin{align*} \mbox{obj. func.: }&V(s_1,s_2,\ldots,s_n) = \sum_{i=1}^{n}{l^T s_i} \rightarrow \min_{x,w,s} \\ \mbox{s.t.: } & \forall i=1,2,\ldots,n: \\ & \qquad s_i \leq s_{\max} x_i \in \R^N \\ & \qquad \forall k=1,2,\ldots,n_l: \\ & \qquad \qquad C w_i^{(k)} = f^{(k)} \in \R^m \\ & \qquad \qquad \sigma_- \bar{s}_i \leq w_i^{(k)} \leq \sigma_+ \bar{s}_i \\ & \sum_{i=1}^{n}{x_i} \leq 1 \\ & x \in \{0;1\}^{nN} ; \qquad 0 \leq s \\ \mbox{with: } & \bar{x}_i = \sum_{j=1, j \ne i}^{n}{x_j} ; \qquad \bar{s}_i = \sum_{j=1, j \ne i}^{n}{s_j} \end{align*} Trivially :math:`0 \leq \sum_{i=1}^{n}{\bar{x}_i} \leq 1` holds and therefore :math:`0 \leq \bar{x}_i \leq 1`. From :math:`x \in \{0;1\}^{nN}` follows :math:`0 \leq \sum_{i=1}^{n}{x_i}` and :math:`s \leq s_{\max}`. The damaged structures are also well defined: .. math:: \bar{s}_i = \sum_{j=1, j \ne i}^{n}{s_j} \leq s_{\max} \sum_{j=1, j \ne i}^{n}{x_j} = s_{\max} \bar{x}_i \leq s_{\max} Clearly :math:`0 \leq \bar{s}_i` holds. So we get for the damaged structures :math:`0 \leq \bar{s}_i \leq s_{\max}` and for the complete structure :math:`0 \leq s \leq s_{\max}`. This is done by :py:func:`ttrr.ttrr_redundant` for redundancy = :math:`n \geq 2`. .. [#] RAID is an acronym for Redundant Array of Independent Disks. http://en.wikipedia.org/wiki/RAID .. [#2] MIP: mixed integer problem; MILP: mixed integer linear program. Because MIP in general is too time consuming everybody is only thinking of MILP and often calls it MIP. Robust Redundant Topology Optimization of Truss ----------------------------------------------- :index:`Robust redundant topology optimization` is nothing more than redundant topology optimization with appropriate load cases. This concept of redundant robust topology optimization is described in [MohrSteinMatziesKnapek2013]_. .. digraph:: succession rankdir=LR; "robust" [label="ttrr.ttrr_robust"]; "perturbaton" [label="if necessary:\ncreate perturbed forces"]; "ttrr" [label="ttrr.ttrr"]; "redundant" [label="ttrr.ttrr_redundant"]; "robust" -> "perturbaton" -> "ttrr"; "perturbaton" -> "redundant"; Uniqueness of the solution -------------------------- .. index:: linear problem; uniqueness, uniqueness of a linear problem For a linear problem possibly more than one edge of the polyhedron exists as a solution. So the solution of a linear problem is not unique --- but the optimal object value is. Let us consider the small example: .. math:: :nowrap: \begin{align*} \mbox{obj. func.: }& z^T x := \left(\begin{array}{cc}1&1\\\end{array}\right) x \rightarrow \max_{x} \\ \mbox{s.t.: } & A x := \left(\begin{array}{cc}1&1\\\end{array}\right) x \leq b := 1 \\ & x \geq 0 \end{align*} A visualization is given in Fig. :ref:`fig-simpleLP`. .. _fig-simpleLP: .. figure:: mathematical_background_simple_lp.* :width: 7cm :alt: visualization of the LP :align: center Visualization of the simple LP In the visualization we see a level curve for :math:`z^T x = 0` and the polyhedron :math:`\{x \in \R^2 : A x \leq b \wedge x \geq 0\}` bounded by the constraints. Without a lengthy calculation we can see that the optimal value of the objective function :math:`z^T x` is :math:`1`. The edges :math:`(0,1)` and :math:`(1,0)` are feasible, optimal and have both the same value of the objective function --- the same for all points between these edges. So in general you cannot expect a unique solution for a linear program. The simplex algorithm will find an edge as a solution. An interior point algorithm will typically find a solution between edges which are optimal. The following two pictures from :doc:`example_2_load_cases_3_dim` show this behavior --- both are optimal for the same problem. +---------------------------------------------------------------------+------------------------------------------------------------------------------+ | simplex | interior point | +---------------------------------------------------------------------+------------------------------------------------------------------------------+ | .. image:: ../ttrr_src/examples_results/2_load_cases_3_dim_apng.png | .. image:: ../ttrr_src/examples_results/2_load_cases_3_dim_interior_apng.png | | :width: 7cm | :width: 7cm | +---------------------------------------------------------------------+------------------------------------------------------------------------------+ For the redundant topology optimization we have in every case more than one solution. For example let :math:`\tilde{x},\tilde{w},\tilde{s}` with :math:`\tilde{s} = (\tilde{s}_1,\tilde{s}_2,\ldots,\tilde{s}_n)` be a solution, then :math:`(\tilde{s}_2,\tilde{s}_1,\ldots,\tilde{s}_n)` is one, too. For the example :doc:`example_2_load_cases_3_dim` for a "1/3"-redundancy this is shown in the following table: +-----------------------------------------------------------------------------+-----------------------------------------------------------------------------+-----------------------------------------------------------------------------+ | .. image:: batch_2_load_cases_3_dim/2_load_cases_3_dim_redundant_1_apng.png | .. image:: batch_2_load_cases_3_dim/2_load_cases_3_dim_redundant_4_apng.png | .. image:: batch_2_load_cases_3_dim/2_load_cases_3_dim_redundant_2_apng.png | | :width: 4.5cm | :width: 4.5cm | :width: 4.5cm | +-----------------------------------------------------------------------------+-----------------------------------------------------------------------------+-----------------------------------------------------------------------------+ | .. image:: batch_2_load_cases_3_dim/2_load_cases_3_dim_redundant_3_apng.png | .. image:: batch_2_load_cases_3_dim/2_load_cases_3_dim_redundant_5_apng.png | .. image:: batch_2_load_cases_3_dim/2_load_cases_3_dim_redundant_6_apng.png | | :width: 4.5cm | :width: 4.5cm | :width: 4.5cm | +-----------------------------------------------------------------------------+-----------------------------------------------------------------------------+-----------------------------------------------------------------------------+ To mostly overcome this behavior we can require an order of the volumes of the structures: .. math:: l^T s_1 \leq l^T s_2 \leq \ldots \leq l^T s_n This is implemented in :py:func:`ttrr.ttrr_redundant` with the flag "sort". Scaling Factors --------------- .. index:: scaling factor To look at scaling behavior of the above linear programs let us consider the following structure: .. math:: :nowrap: \begin{align*} \mbox{obj. func.: }&V(s) = l^T s \rightarrow \min_{s,w} \\ & C w = f \\ & \sigma_- s \leq w \leq \sigma_+ s \\ & 0 \leq s \leq s_{\max} \end{align*} Material, Load and Design Space ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ The :index:`material parameters` for the topology optimization of truss given to :py:func:`ttrr.ttrr`, :py:func:`ttrr.ttrr_robust` or :py:func:`ttrr.ttrr_redundant` are: * :math:`s_{\max}` : the maximal allowed cross section of a bar * :math:`\sigma_-` : lower limit for the allowed stress, e.g. compressive strength * :math:`\sigma_+` : upper limit for the allowed stress, e.g. tensile strength * density : density The density is only used to calculate the mass of the structure out of its volume. :math:`s_{\max}` is necessary for :py:func:`ttrr.ttrr_redundant` and gives an additional feature to :py:func:`ttrr.ttrr` and :py:func:`ttrr.ttrr_robust`. If its not wanted a huge value would disable it. But in any case it gives with :math:`0 \leq s \leq s_{\max}` bounding values. So we have a closed feasible region. :math:`\sigma_-` and :math:`\sigma_+` are really used material parameters. Typical topology optimization is a preliminary design. Therefore it must not be totally accurate to the praxis. Typically :math:`\sigma_- = -\sigma_+ < 0` is chosen. So we can assume :math:`\sigma_- = -\sigma_+ = -\sigma` for some :math:`\sigma > 0`. .. math:: :nowrap: \begin{align*} \mbox{obj. func.: }&V(s) = l^T s \rightarrow \min_{s,w} \\ & C w = f \\ & -\sigma s \leq w \leq \sigma s \\ & 0 \leq s \leq s_{\max} \end{align*} Easily to see the following program results in an equivalent optimum. With :math:`\tilde{s} = \sigma s` we get: .. math:: :nowrap: \begin{align*} \mbox{obj. func.: }&V(\tilde{s}) = l^T \tilde{s} \rightarrow \min_{\tilde{s},w} \\ & C w = f \\ & -\tilde{s} \leq w \leq \tilde{s} \\ & 0 \leq \tilde{s} \leq \frac{1}{\sigma} s_{\max} \end{align*} This is only a scaling. The magnitude of the loads is only a scaling, too. With :math:`\tilde{w} = \alpha w` for some :math:`\alpha > 0` we get: .. math:: :nowrap: \begin{align*} \mbox{obj. func.: }&V(\tilde{s}) = l^T \tilde{s} \rightarrow \min_{\tilde{s},\tilde{w}} \\ & C \tilde{w} = \alpha f \\ & -\alpha \tilde{s} \leq \tilde{w} \leq \alpha \tilde{s} \\ & 0 \leq \tilde{s} \leq \frac{1}{\sigma} s_{\max} \end{align*} And furthermore with :math:`\hat{s} = \alpha \tilde{s} = \alpha \sigma s` the following program is equivalent, too: .. math:: :nowrap: \begin{align*} \mbox{obj. func.: }&V(\hat{s}) = l^T \hat{s} \rightarrow \min_{\hat{s},\tilde{w}} \\ & C \tilde{w} = \alpha f \\ & - \hat{s} \leq \tilde{w} \leq \alpha \hat{s} \\ & 0 \leq \hat{s} \leq \frac{1}{\alpha \sigma} s_{\max} \end{align*} Keep in mind it is not required to change the objective function. For :math:`\beta > 0` the minimum of the functions :math:`g(x)` and :math:`\beta g(x)` is reached at the same point. The size of the design space is not found directly in the linear program. In the geometry matrix :math:`C` occur only the angles of the bars. So the magnitude of the size of the design space in meter or millimeter will result in the same optimal structure, because :math:`l` is only in the objective function and will be changed by a scalar factor. Numeric ^^^^^^^ A scaling of .. math:: :nowrap: \begin{align*} \mbox{obj. func.: }&V(s) = l^T s \rightarrow \min_{s,w} \\ & C w = f \\ & \sigma_- s \leq w \leq \sigma_+ s \\ & 0 \leq s \leq s_{\max} \end{align*} do not change the (relative) condition of the problem. But it can change the absolute condition and therefore helps the computer by the calculation. To do so the functions :py:func:`ttrr.ttrr`, :py:func:`ttrr.ttrr_robust` or :py:func:`ttrr.ttrr_redundant` have the parameters :math:`\alpha =` scale_sigma and :math:`\beta =` scale_forces. To get the picture let us see how they transform the above problem: .. math:: :nowrap: \begin{align*} \mbox{obj. func.: }&V(s) = l^T s \rightarrow \min_{s,w} \\ & C \beta w = \beta f \\ & \alpha \sigma_- \frac{\beta}{\alpha} s \leq \beta w \leq \alpha \sigma_+ \frac{\beta}{\alpha} s \\ & 0 \leq s \leq s_{\max} \end{align*} Now the calculation will be done with :math:`\tilde{s} = \frac{\beta}{\alpha} s`, :math:`\tilde{w} = \beta w`, :math:`\tilde{f} = \beta f`, :math:`\tilde{\sigma}_- = \alpha \sigma_-`, :math:`\tilde{\sigma}_+ = \alpha \sigma_+` and :math:`\tilde{s}_{\max} = \frac{\beta}{\alpha} s_{\max}`: .. math:: :nowrap: \begin{align*} \mbox{obj. func.: }&V(\tilde{s}) = l^T \tilde{s} \rightarrow \min_{\tilde{s},\tilde{w}} \\ & C \tilde{w} = \tilde{f} \\ & \tilde{\sigma}_- \tilde{s} \leq \tilde{w} \leq \tilde{\sigma}_+ \tilde{s} \\ & 0 \leq \tilde{s} \leq \tilde{s}_{\max} \end{align*} Afterwards :math:`\tilde{s}` and :math:`\tilde{w}` will be transformed back. For a maximum magnitude of :math:`1` of an element from :math:`\tilde{f}` the absolute condition of the equation will be a little bit better. For a maximum magnitude of :math:`1` of an element from :math:`\tilde{\sigma}_-` and from :math:`\tilde{\sigma}_+` the condition of the inequation will be significant better. If :math:`\alpha = \beta` the design variables :math:`s` and :math:`s_{\max}` resp. will not be changed. .. only:: html References: .. [Ben-TalElGhaouiNemirovski2009] Aharon Ben-Tal ; Laurent El Ghaoui ; Arkadi Nemirovski: Robust Optimization. Princeton University Press, 2009. http://press.princeton.edu/titles/9099.html .. [BendsoeSigmund2004] Bendsoe, Martin P. ; Sigmund, Ole: Topology Optimization. Berlin, Heidelberg : Springer, 2004. ISBN 3-540-42992-1 .. [Marti2003] Kurt Marti: Stochastic optimization methods in optimal engineering design under stochastic uncertainty. In: ZAMM, 83 (12): 795--811, 2003. http://dx.doi.org/10.1002/zamm.200310072 .. [Michell1904] Michell, Anthony George M.: The limits of economy of material in frame structures. In: Phil. Mag. 8 (1904), Nr. 47, S. 589-597. http://dx.doi.org/10.1080/14786440409463229. .. [Mohr2011] Mohr, Daniel P.: Redundante Topologieoptimierung. Neubiberg, Universitaet der Bundeswehr Muenchen, Fakultaet fuer Luft- und Raumfahrttechnik, Diss., Dezember 2011. http://nbn-resolving.de/urn/resolver.pl?urn=urn:nbn:de:bvb:706-2664 .. [MohrSteinMatziesKnapek2012] Mohr, Daniel P. ; Stein, Ina ; Matzies, Thomas ; Knapek, Christina A.: Robust Topology Optimization of Truss with regard to Volume. In: arXiv - Mathematics, Optimization and Control (2012). http://arxiv.org/abs/1109.3782v2 .. [MohrSteinMatziesKnapek2013] Mohr, Daniel P. ; Stein, Ina ; Matzies, Thomas ; Knapek, Christina A.: Redundant robust topology optimization of truss. In: Optimization and Engineering (2013), 1-28. http://dx.doi.org/10.1007/s11081-013-9241-7 .. [PattersonGibsonKatz1987] Patterson, David A. ; Gibson, Garth A. ; Katz, Randy H.: A Case for Redundant Arrays of Inexpensive Disks (RAID). Version: Dec 1987. 1987 (UCB/CSD-87-391). --- Forschungsbericht. --- EECS Department, University of California, Berkeley. http://www.eecs.berkeley.edu/Pubs/TechRpts/1987/5853.html .. [Przemieniecki1968] Janusz S. Przemieniecki: Theory of matrix structural analysis. McGraw-Hill, Inc., New York, 1968. .. [Shannon1948] Shannon, Claude E.: A Mathematical Theory of Communication. Version: 1948. 1948 (1). --- Forschungsbericht. --- 379-423, 623-656 S. --- The Bell System Technical Journal. http://cm.bell-labs.com/cm/ms/what/shannonday/paper.html