Game theory approach to robust topology optimization with uncertain loading

The paper concerns robustness with respect to uncertain loading in topology optimization problems. Using a game theoretic framework we formulate problems, or games, defining generalized Nash equilibria. In each game a set of topology design variables aim to find an optimal topology, while a set of load variables aim to find the worst possible load. Several numerical examples with uncertain loading are solved in 2D and 3D. The games are formulated using global stress, mass and compliance as objective functions or constraints.


Introduction
We consider structural topology optimization (TO) problems where there is uncertainty regarding the loads that the optimized design will be subjected to. Load uncertainty is typically the case for parts in a fighter aircraft, which, due Erik Holmberg erik.holmberg@liu to advanced maneuvers, may have high accelerations in any direction. It is well known that an optimized structure may have very bad performance if subjected to a load which it was not optimized for; therefore, even if we would use a very large number of load cases in the optimization, a small perturbation of one of those loads may be catastrophic for the structural integrity.
In this work, rather than using a limited number of fixed load cases, the load may be any load in an infinite uncertainty set defined by known maximum loads and assuming an elliptical variation in space. This type of problem falls into the field of robust optimization (RO) (Ben-Tal et al. 2009), where one often distinguishes between stochastic and deterministic approaches. In the case of TO the former includes reliability-based design optimization (RBDO), where so-called reliability indices (Valdebenito and Schuëller 2010), quantifying the probability of structural failure, are used as objective or constraints, and robust design optimization (RDO) methods, where one minimizes expected values and standard deviations of, e.g., compliance given loads from some probability distribution (Evgrafov et al. 2003;Dunning and Kim 2013;Jansen et al. 2015). In deterministic, or worst-case, RO one assumes nothing about the probability distribution of, in this case, the load-variation, except upper and lower bounds on e.g. magnitudes, and this is what is being done in the present paper.
Deterministic, robust TO with uncertain loading has previously been studied for minimum compliance formulations. These problems may, assuming small deformations and using certain load parametrizations, be cast as generalized eigenvalue problems (Brittain et al. 2012;Cherkaev and Cherkaev 2008;Takezawa et al. 2011) or as semidefinite programming problems ; Thore et al. 2015). In this paper however we propose a much more general game theoretic framework for TO under loaduncertainty including a wide range of different objective functions and constraints.
Very briefly, game theory defines a problem class where two (or more) players, given certain strategies (or variables), cooperate or are in conflict while trying to optimize their profit (Aubin 1979). Game theory approaches for structural optimization have been used in a number of papers, but only two papers (known by the authors) involve TO: Habbal (2005) formulated a TO framework to model growth of cancerous tumors and Habbal et al. (2004) optimized a thermoelastic system where a heat source as well as an external load were applied to the design space; one player was aiming to minimize the compliance of the thermoelastic system by varying the topology, the other player was aiming to minimize the temperature in the structure by adding structural elements, such as fins, in order to maximize the heat flow to the surroundings.
References to other fields of structural optimization where game theory has been applied include Kobelev (1993) where size optimization was performed on truss structures with varying loads. Using the same convex pay-off functional for both players, it was shown that the solution to the game was obtained as an eigenvalue problem. Périaux et al. (2001) used genetic algorithms in a shape optimization problem to find a Nash equilibrium for players optimizing the shape and flow in a nozzle. Banichuk (1973) calculated analytical solutions to optimization problems of elastic beams subjected to loads of various forms from a predefined set. He noted that for some classes of problems there exists a unique worst-case load, for which it is possible to find an optimal design that is optimal for the entire set, and also for that unique load. But for some problems the worst-case loading is not unique; consequently, an optimized design will not be optimal for any of the loads in the set if these are applied as unique load cases. This observation is related to the multiple eigenvalue problem discussed in Holmberg et al. (2015). Uncertainty, not only concerning loading, but also elastic moduli and material defects such as cracks was considered by Banichuk and Neittaanmäki (2007).
For the special case when both players use the same objective function g, known as a zero-sum game, Aubin and Ekeland (1984, Proposition 1, Chapter 6) show that a non-cooperative equilibrium is also a solution to a min-max problem with g as objective (cf. Appendix A). Choosing the compliance as the objective function and a suitable loading parametrization one can retrieve the generalized eigenvalue problems (Brittain et al. 2012;Cherkaev and Cherkaev 2008;Takezawa et al. 2011;Kobelev 1993) or semi-definite programming problems Thore et al. 2015) mentioned above. When applicable these formulations may be more efficient, but compared to the proposed game theory framework they are very limited in that they only apply to zero-sum games with certain choices of functions and parametrizations.
In order to find designs that are robust, and optimized, with respect to an uncertain load we here seek so-called generalized Nash equilibria (GNE), 1 which defines situations where none of the players have an incentive to change strategy, unless the other player does so. In our games we parametrize the load vector using one set of variables that control the load within the uncertainty set. We also have the standard TO variables that determine which elements in the finite element-discretized design domain that represent material and which represent holes. Both sets of variables influence the state equation and are chosen such that they are in conflict with each other; in game theory vocabulary we say that we have a two-player non-cooperative game (Aubin 1979).
The variables and the linear elastic model as well as the game theoretic formulation in generic form is introduced in Section 2. Optimality conditions are then given in Section 3. The algorithm that we use is stated in Section 4, and the design and load parametrization is described in Section 5. Three specific functions (to be used as objective or constraint functions) are given in Section 6. In Section 7 we discuss the choice of additional load cases in our games to obtain convergence by the proposed algorithm. Numerical solutions to three concrete instances, in both 2D and 3D, of the generic game are presented in Section 8.

Structural model in the game theoretic framework
We formulate games with two sets of design variables and two objective functions. The first set of variables is the topology variables x and the second set is the load variables θ . The topology variables x are chosen to find the best possible structure by minimizing the function g 1 (x, θ ), which in this paper is stiffness, mass or stress. The load variables θ are chosen to find the worst possible loading by maximizing the function g 2 (x, θ ), for which there are also several possible choices; we use stiffness and stress in this paper. The topology variables are restricted to where c e defines upper bounds on values of the k constraint functions c e . The m design variables are subject to box constraints where 1 implies that the corresponding element contains material and the small number ε > 0 that it is empty. The load variables are restricted to in which s is the number of spatial dimensions. The box constraints in can be used to restrict the load variation, but in this paper they are chosen sufficiently loose so that they do not become active and are included only to improve numerical performance.
Both sets of variables influence the state equation where K(x) ∈ R n×n is the global stiffness matrix, and U = u 1 , . . . , u f ∈ R n×f contains the nodal displacements obtained for the f load cases in

The game in generic form
Following the terminology of game theory (Aubin 1979), a game is characterized by players who have certain objective functions (sometimes called e.g. loss-or pay-off functions) which they influence using strategies. Our problem is a twoplayer non-cooperative game, non-cooperative meaning that the two players can change only their own strategy. The first player wants to minimize g 1 (x, θ ), using as strategy the design variables x. The second player wants to maximize g 2 (x, θ ) and uses the load variables θ as strategy. The sets X (θ) and introduced above are now referred to as strategy sets of the respective player.
The meaning of the term solution to a game can be ambiguous. The most commonly used solution concepts are Pareto optima and (generalized) Nash (or non-cooperative) equilibria. The games studied in this paper are all of zerosum character, 2 and for a zero-sum game, every point is in fact a Pareto optimum, making this solution concept unsuitable. A solution to our game is thus a generalized Nash equilibrium, i.e. a point (x * , θ * ) satisfying Problems (2a) and (2b) are in general non-convex (resp. non-concave), large-scale optimization problems for which computing globally optimal solutions is currently not practical. Therefore one should think here of GNE as meaning, a priori, local GNE, with x * and θ * locally optimal solutions to (2a) and (2b), respectively.
An important question is of course whether the game (2) admits any equilibria at all? Even if (2) was a standard and not a generalized Nash game, the classical Nash existence theorem (Aubin 1979, p. 267) would not be applicable since, e.g., g 1 is not, in general, convex and g 2 not concave. However, absent a mathematical proof, the numerical results shown below provide strong evidence of existence of solutions for our particular examples, and we note that this difficult topic has been the subject of recent research (Pang and Gesualdo 2011).

First order necessary optimality conditions
We introduce two Lagrangian functions where λ e , ξ i , η i , κ i and γ i are Lagrange multipliers.

Algorithm for finding a generalized Nash equilibrium
An overview of algorithms for finding GNEs is given by Facchinei and Kanzow (2010). Unfortunately such algorithms are much less developed, in terms of both theory and implementation, than, e.g., those for non-linear optimization problems, for which several high-quality codes, often globally convergent under reasonable assumptions, are available.
One way to find a GNE is to solve the optimality system (3a)-(3q) directly using Newton-based methods (Christensen et al. 1998;Dreves et al. 2011;Facchinei and Kanzow 2010). Such methods, however, require second derivatives, or at least approximations thereof. In particular the Hessian of L 1 is large and dense and Newton-based methods are therefore expected to be too expensive, unless perhaps (limited memory) quasi-Newton approximations can be used. Périaux et al. (2001) and Habbal et al. (2004) used a nonlinear Jacobi-type algorithm (Facchinei and Kanzow 2010) where the two players update their strategies simultaneously, i.e. both optimization problems are solved at the same time and independent of each other. This algorithm has the advantage that it allows for computing the two strategies in parallel. However, in our game the two players' strategies are highly counteracting in (1) and we have found it more efficient to update the strategies in sequence. Therefore we suggest the following nonlinear Gauss-Seidel-type (Facchinei and Kanzow 2010) algorithm: In the examples in Section 8, the problems in Step 1 and Step 2 are considered solved whenever the change in the objective function value is below a given threshold and the constraints are satisfied within a certain tolerance. As for the stopping criteria in Step 3, Algorithm 1 is terminated whenever only one iteration of the TO problem in Step 2 was required in three successive outer iterations (an outer iteration consists of Step 1, 2 and 3).
As pointed out by Facchinei and Kanzow (2010), if Algorithm 1 converges, it does so to a GNE (it suffices that all functions involved are continuous). The conditions under which convergence occurs is however not clear. The most promising approach towards ensuring convergence for a wide class of problems seems to be to include so-called proximal terms in the objective of the respective problem to penalize too rapid changes of the strategies (Attouch et al. 2008;Facchinei and Kanzow 2010). This has not been done here; instead we propose in Section 7 to modify the games themselves such that convergence is obtained when we apply Algorithm 1.

Topology variables
We use a standard SIMP-approach (Bendsøe 1989) for the TO. The topology design variables x are filtered using a linear design variable filter (Bruns and Tortorelli 2001) that may be written in which, for element j , e j is its centroid and v j its volume; the filter radius is denoted R and || · || is the Euclidean norm. We denote ρ i (x), i = 1, . . . , m,as physical variables as they are used to define the properties of the structure. Using the SIMP-approach, the stiffness matrix reads where K i is an expanded element stiffness matrix and p > 1 is a penalization exponent. The lower bound ε introduced in X (θ) ensures that the stiffness matrix is positive definite for every feasible design x, assuming the structure is appropriately supported. The filter in (4) results in a transition zone between material and holes (Sigmund 2007). Due to the penalization in (5) the intermediate design variable values in this transition give a non-physical stiffness, which makes them undesirable. Therefore, in order to obtain a more "black-and-white" design without the intermediate design variable values, we remove the filter after convergence of Algorithm 1 and continue with a TO using only those ρ i (x) that are close to a boundary as design variables, while those that are at the  Holmberg et al. (2015) for additional details). As the design changes are typically quite small we have chosen not to update the loads after the filter has been removed, and we have found (3k)-(3q) to be very close to satisfied anyway, see Figs. 6b, 9c and 12c. However, we note that there is no guarantee that the final design is a GNE. In that sense, we consider the design obtained after Algorithm 1 to be the optimized design and the final design the result of a post-processing step that, at least, satisfy (3a)-(3j).

Load variables
We consider loads varying in an uncertainty set where the maximum loads in all directions are known, but the direction of the load may vary. Such loads can for example occur due to accelerations in an aircraft, where there are restrictions on the allowable accelerations in different directions, but accelerations may occur in any direction depending on the maneuver. The loading, for load case , reads where the uncertainty vector r ∈ R s is a unit vector and L ∈ R s×n a matrix containing the maximum loads. This type of loading has been used in Holmberg et al. (2015) where the components of r were used to vary the loading. In this paper we choose to parameterize r by the angles θ = (θ, φ) T in a spherical coordinate system (Fig. 1), reducing to a polar coordinate system for 2D-problems (φ = π/2). These angles are used as variables to define the direction of the load.
The first load case f 1 (θ), i.e. the first column of F (θ) in (1), is for a 3D-problem expressed as The matrix L is built by a number of s × s diagonal blocks, but only one such block is written explicitly in (6). Each block corresponds to the degrees of freedom of one node; if no load is applied at that node the block contains zeros, otherwise the maximum loads x max , y max and z max are placed on the diagonal. For notational simplicity we have assumed that the maximum loads are given in the global coordinate directions and we assume that the maximum loads are the same in all nodes, i.e. no further index is added.
Note that we expect the worst load to occur on the ellipsoid (or ellipse in 2D), not for a load of lower magnitude inside it. For the functions introduced in Section 6 this is proved mathematically using the fact that the maximum value of a convex function is found at the extreme points of the feasible set (Rockafeller 1972, Corollary 32.3.1). The Fig. 2 Load variable and load for a 2D-problem convexity of both the compliance and the stress measure introduced in Section 6 as functions of the load follow by viewing each as composed of a non-decreasing and a convex function (Boyd and Vanderberghe 2004, p. 97).

Extensions to more general loading scenarios
It is possible to extend (6) to more general loading scenarios: -If there are several external loads varying with the same θ it is straightforward to add maximum loads in the respective block of L and each block may have different maximum loads on the diagonal. -If there are several load cases, with possibly different maximum loads or loads applied at different positions, we add another columnL T r (θ) to F (θ), whereL is built in the same way as L and the number of variables in θ is increased accordingly. -Additional load cases that vary with the same θ and are restricted by the same L are introduced by adding columns in r(θ). This is done in this paper, as described in Section 7.

Objective functions, constraints and gradients
The game theoretic formulation allows for a wide range of choices regarding g 1 , g 2 and c e in (2). Since we do not have an existence proof we cannot be too specific, but the possible choices is at least restricted by the requirement that each of the problems in (2) be well-posed. For this it suffices (Andréasson et al. 2005, Theorem 4.7) that, e.g., X (θ) and be non-empty and compact and g 1 and g 2 lower, respectively upper, semi-continuous, conditions that are not difficult to see to in practise. For numerical efficiency it is preferable that g 1 , g 2 and c e be at least differentiable so that gradient-based optimization algorithms can be applied.
In this paper we consider three functions: compliance, global stress, and mass, which we define and differentiate below.

Compliance
The compliance for load case is where u (x, θ ) is part of the solution to (1).
The gradient of (7) with respect to x can be found in e.g. Christensen and Klarbring (2008). The derivative with respect to θ t , t = 1, . . . , (s − 1) is given by where the second and third step follow from (1). From (6) we find that where ∂r (θ)/∂θ t is the derivative of the simple trigonometric functions in (6). The θ t -derivative of (7) therefore finally reads

Stress
A global stress measure, σ G (θ), = 1, . . . , f , is created for each load case using a P -norm, with exponent P ≥ 1, of the vector of local von Mises stresses (Holmberg et al. 2013b): where σ vM a (x, θ ) is the von Mises stress for load case in stress evaluation point a and d is the number of such points. The von Mises stresses are based on penalized stresses and the penalized stress vector (in Voigt notation) in stress evaluation point a for load case is calculated as where E is the constitutive matrix, B a is the expanded strain-displacement matrix corresponding to stress evaluation point a and index i denotes the element that this point belongs to. The stress penalization ρ i (x) q , with 0 < q < 1, was suggested by Bruggi (2008) and also used successfully by e.g. Le et al. (2010) and Holmberg et al. (2013a). The gradient of the global stress measure (9) with respect to x is given in Holmberg et al. (2013bHolmberg et al. ( , 2013a and is not repeated here. The derivative of (9) with respect to θ t reads The first factor in (11) is obtained from (9) as the second factor follows trivially from the definition of the von Mises stress; and the third factor follows from (10) and reads where ∂f (θ)/∂θ t is given in (8). Rather than forming the inverse of the stiffness matrix explicitly we compute K −1 (x)∂f (θ)/∂θ t by solving an additional linear system, at the same cost as solving (1).

Mass
The mass M(x) is where M i is the element mass. The derivative with respect to design variable x b reads where ib was defined in (4). Obviously, since the mass is not a function of the load, it does not make sense to use mass for g 2 in (2b).
Remark 1 Using the global stress measure (9) in a stress constraint we are guaranteed that all local stresses (in the evaluation points in the FE model) are lower than or equal to the stress limit (Holmberg et al. (2013a(Holmberg et al. ( , 2013b) and we consider the stresses in the final designs in Figs. 12b and 14b to be sufficiently close to the stress limit given that TO is essentially a conceptual design tool. If more control over the local stresses is desired, an alternative is to use clustered stress measures where the stress evaluation points are arranged into a number of clusters and one stress measure is calculated for each such cluster (Holmberg et al. 2013a). The computational cost for this is of course higher than for the global approach used here.

Adding load cases for stability
A structure optimized for only one load will typically have bad performance for other directions of the loading. This can cause difficulties in the form of oscillations between Step 1 and 2 when Algorithm 1 is applied to games such as those studied here. This issue is illustrated next by an example, which also suggest how to modify problem (2a) to get convergence for our games.
If the design space is sufficiently large the worst load found at iteration k + 1 will be a load where the defining uncertainty vector r is orthogonal to that found at iteration k; this is shown in Appendix A for a problem where both g 1 and g 2 is compliance. This means that the topology x k+1 can be totally different from x k and this causes an oscillation where the Fig. 3 An example of an intermediate design step, just before the load is about to be updated. a: Only one load case and obviously not a robust design, b: Two load cases (shown simultaneously) and a more robust design load varies and the topology oscillates between different designs. An example where both g 1 and g 2 is compliance is visualized in Fig. 3a, which shows an intermediate (not converged) design where the load direction is just about to be updated. Here the new worst load direction will be orthogonal to the current and major topological changes are required in order to minimize the compliance for the new load.
One way to avoid these oscillations is to add load cases to the TO problem (2a). Motivated by Appendix A, loads defined by orthogonal uncertainty vectors should be a good choice for additional loads. For the 3D-examples in Section 8 however, we found that Algorithm 1 converged faster for orthogonal loads than for orthogonal uncertainty vectors, which is why orthogonal loads are added in this paper. If x max = y max (= z max in 3D) the same loads are obtained by using orthogonal uncertainty vectors as orthogonal loads; but if the loading region is elliptical, orthogonal loads will span the ellipsoid better than the loads due to orthogonal uncertainty vectors. Figure 3b shows a design, in an intermediate design step, obtained using two load cases: the primary load f 1 and a second load f 2 , orthogonal to f 1 .
Adding more load cases might make the algorithm converge faster, but in each iteration we need to solve for more right hand sides in (1) and, if applicable, in the adjoint gradient calculation. From numerical tests we have found it efficient to add loads obtained by rotating (the 2 or 3 orthogonal loads) π/4 radians in each plane; see Fig. 4 for a visualization in 2D. This means that we use four load cases in 2D problems and nine load cases in 3D problems and all loads are defined by θ . The additional computational cost due to the extra right hand sides is usually compensated for by faster convergence. One way of calculating the additional loads is given in Appendix B.
If the objective function in the topology problem (2a), depends on θ we use the sum of the objective values for the load cases and constraints are added if a constraint function Fig. 4 Four load cases used in a 2D problem depends on θ . In the load problem (2b), we use only the first load case (6) -in this way we find the worst load rather than the worst weighted loads.

Examples
There are several possible games that fit into the generic formulation (2). Based on the functions defined in Section 6 we formulate three games that we consider numerically, using two different geometries.
The first geometry is the two-dimensional L-beam seen in Fig. 5a. The L-beam is a challenging test example because of the internal corner with a stress singularity; the stress is very high initially, so the initial design is very far from feasible when we have a stress constraint. The L-beam has outer dimensions 200 × 200 mm, thickness 1 mm and is rigidly attached at the upper boundary. The applied load is due to acceleration of a 12 kg mass attached in the right corner. The maximum accelerations are 10 times the acceleration of gravity: giving a maximum load of 1177 N in any direction in the xy-plane, i.e. x max = y max = 1177  (6). The design domain is discretized with 6400, equal sized, 8-node quadrilateral elements and has 39040 degrees of freedom.
The second geometry is an attachment bracket modelled in 3D, Fig. 5b, with outer dimensions 230 × 60 × 68 mm. The bracket is rigidly attached in one end and attaches to some equipment that weighs 22 kg in the lower part of the other end. The bracket shall be dimensioned for 20 times the acceleration of gravity in the z-direction and 10 % of that in other directions for robustness. The maximum loads are thus x max = y max = 431.6 and z max = 4316 and the load is modelled as a point load applied in the centre of gravity of the equipment, from which stiff rods connect to four attachment points in the bracket. The bracket is modelled with linear eight node hexahedron-and six node pentahedron elements of varying size; in total 25864 elements and 86187 degrees of freedom. The filter radius is set to To avoid numerical issues due to the stress concentration at the point where the load is applied, the elements in the vicinity of that point in the 2D example and in the vicinity of the rods in the 3D example are non-design elements for which the stress is not part of the global stress measure.
The topologies are plotted in a gray-scale where white implies void, black solid material and gray is elements with intermediate design variable values. The 3D topologies have gray element lines for ease of visualization. The stress plots should be viewed in color. The design material is an aluminium with Young's modulus 71000 MPa, Poisson's ratio 0.33 and density 2.8 × 10 −9 tonne/mm 3 .
In all numerical examples the initial design is an equal distribution of material; the lower bound on the design variables is set to ε = 0.001; the SIMP and stress penalization exponents to p = 3 and q = 0.5, respectively; and the Pnorm exponent to P = 24. The local penalized stress (10) is evaluated in the element centroid.
The optimization problems in Step 1 and 2 of Algorithm 1 are solved using MMA, the Method of Moving Asymptotes (Svanberg 1987). All examples have converged at a low number of outer iterations k in Algorithm 1, typically in the order of 5 to 15, and the number of inner iterations in Step 1 and 2 is drastically decreased with increasing k, the last iterations only involving one inner iteration. The total number of inner iterations has been between approximately 200 and 1500, using conservative settings of the allowable updates in MMA and fine convergence criteria. As the load problem in Step 1 has only s − 1 variables and simple box constraints, it is solved in a few iterations. Thus, the main computational cost is due to solving the TO problem in Step 2.

Game formulation 1: compliance
In the first formulation we minimize compliance, subjected to a mass constraint, using the topology variables and search for the load giving the highest compliance using the load variables, i.e.
The design seen in Fig. 6  As a reference, the design obtained for non-robust optimization with a fixed load, f 1 (θ = 1.5π ), is given in Fig. 7. Comparing the robust and non-robust designs we find that the main difference is the diagonal structural member between the two vertical members, which increases the stiffness for a load in the x-direction.
The curve in Fig. 6b shows g 2 (x * , θ ) = C 1 (x * , θ ) for θ ∈ [0, 2π ] and the dashed vertical line shows θ = θ * . We find that f 1 (θ * ) is actually the load that gives the maximum compliance, showing that the design in Fig. 6a is a robust design.
The 3D design obtained using (G 1 ) and the upper bound M = 0.1M is seen in Fig. 8.

Game formulation 2: stress minimization
In the second formulation we minimize the weighted global stress, subjected to a mass constraint, using the topology variables and maximize the global stress using the load variables: where X 2 = X 1 . For the allowable mass M = 0.45M we obtain the design in Fig. 9a. A smooth radius has been created in the internal corner in order to avoid the stress singularity and the thickness of the structural parts are dimensioned such that the stress in the structure is minimized considering all possible loads; therefore, some parts will have low stress for (c)  Fig. 9b shows the von Mises stresses for the load f 1 (θ * ). Studying the stress plot it is obvious that some structural parts have very low stress, implying that they are not necessary for this particular load case.
Again, the curve in Fig. 9c shows that θ * is, within some tolerance, the maximizer of the global stress measure g 2 (x * , θ ) = σ G 1 (x * , θ ), implying robustness. A non-robust design, using the fixed load f 1 (θ = 1.5π ) and the same M, is seen in Fig. 10. In the non-robust design the structural parts are thicker as they are only dimensioned for one specific load. The stress is thus lower for that load case than what it is for the robust design, but higher for other load cases.
The equilibrium design for the 3D example using M = 0.1M is shown in Fig. 11a and the von Mises stresses for the first load case are shown in Fig. 11b.

Game formulation 3: stress constraint
In (G 2 ) we used a fixed limit on the available mass, which in the 2D example was found manually such that the optimized structure had a maximum von Mises stress approximately equal the yield limit of the material, 350 MPa. In formulation 3 we are instead aiming for the lightest design that satisfy a stress constraint, and using the load variable we seek the direction giving the highest global stress. This formulation reads The equilibrium design obtained using (G 3 ), seen in Fig. 12a, has, as for (G 2 ), a radius in the internal corner so that the stress singularity is avoided, and the structural members are dimensioned such that the global stress is lower than the constraint limit, chosen equal to the yield limit of the material, σ G = 350 MPa. The maximum local stress for  Fig. 12b, is 322 MPa. The reason why the maximum local stress is slightly lower than σ G is because the global stress measure (9) overestimates the maximum local stress (Holmberg et al. 2013b).
Comparing again with a non-robust design, optimized for the fixed load f 1 (θ = 1.5π ) and seen in Fig. 13, we find that the non-robust optimization has created fewer structural (c) members with approximately the same stress for this particular load, whereas the robust design has more structural members and not all of these are fully stressed for the load f 1 (θ * ). However, when varying θ all parts should become fully stressed for some specific load.
The curve in Fig. 12c shows that we are at a local maxima for the load problem, but that the global maxima is only about 2 % higher. As a comparison we also plot the same curve for the non-robust design, Fig. 13c, where we find that the maximum global stress is approximately 450 % higher Mises stress for the first load case than the stress limit for a load orthogonal to that it was optimized for.
For the 3D example optimized using (G 3 ) we find that the optimized design in Fig. 14a is simple and easy to interpret and that the maximum von Mises stresses, Fig. 14b, are below the stress limit σ G = 350 MPa. However, some structural members have become very thin compared to the element size and some gray elements remain. This implies that, for this load and stress constraint, a finer mesh should have been used if better precision was required from the model. Now the stress constraint is satisfied and it is not possible to remove more elements without removing a structural member.

Concluding remarks
A game theory approach to robust TO with uncertain loading has been developed and exemplified using three different games for design of both 2D and 3D structures.
The nature of the proposed non-cooperative games, between the structure and the external loads, is such that convergence is difficult to obtain -an element may be very important for some loads but completely unnecessary for others, and this typically leads to oscillations in the design variable values. By addition of certain load cases (see Section 7) such oscillations can be avoided, thus making Algorithm 1 efficient. The efficiency is demonstrated by the small number of outer iterations and the almost black-andwhite optimized designs. The few gray elements remaining are in part explained by the filtering technique combined with the varying element size in the 3D mesh, but also by the discretization as discussed in connection with Fig. 14.
The proposed game theoretic framework allows formulation of a wide class of relevant structural optimization problems, most of which cannot be cast as semi-definite programming problems or optimization problems with eigenvalues. The numerical examples give an idea of the broad range of problems that can be given in the game theoretic framework; for example that we can apply stress constraints despite the load uncertainty. in order to find the worst loading at iteration k. The