Robust Design of a Smart Structure under Manufacturing Uncertainty via Nonsmooth PDE-Constrained Optimization

We consider the problem of finding the optimal shape of a force-sensing element which is integrated into a tubular structure. The goal is to make the sensor element sensitive to specific forces and insensitive to other forces. The problem is stated as a PDE-constrained minimization program with both nonconvex objective and nonconvex constraints. The optimization problem depends on uncertain parameters, because the manufacturing process of the structures underlies uncertainty, which causes unwanted deviations in the sensory properties. In order to maintain the desired properties of the sensor element even in the presence of uncertainty, we apply a robust optimization method to solve the uncertain program.The objective and constraint functions are continuous but not differentiable with respect to the uncertain parameters, so that existing methods for robust optimization cannot be applied. Therefore, we consider the nonsmooth robust counterpart formulated in terms of the worst-case functions, and show that subgradients can be computed efficiently. We solve the problem with a BFGS--SQP method for nonsmooth problems recently proposed by Curtis, Mitchell and Overton.


Introduction
Mechanical structures such as bridges, cars or airplanes are meant to bear the occurring loads and to maintain their intentional functionality throughout their lifetime. Failures, accidents and breakdowns have to be avoided by all means. However, there is an inherent uncertainty throughout the whole life span of a system or a component. Assumptions made during the design phase may differ from the actual usage scenario of systems and components. Theoretical calculated values may vary from the loads that actually occur, e.g., due to corrosion or aging phenomena. After the design the system has to pass the production phase. Here, additional uncertainty in the manufacturing processes occurs, which leads to uncertain properties of the actual parts due to geometrical deviations or unknown material properties. Consequently, uncertain parts are produced which are going to show an unknown behavior during the usage phase. This situation cannot be accepted for safety critical structures. Thus, many measures are taken throughout each phase, such as strict quality control and real-life tests. One additional method to approach this issue is to identify the loads and states of the structures that actually occur during the usage phase. To mitigate the effects of uncertainty by monitoring the states and loads of parts [14] and processes [12], several methods are known, such as Structural Health Monitoring or Process Monitoring. By identifying forces it is possible to react to unforeseen events or to shut down systems before reaching safety critical conditions. Within the SFB 805, a rod framework is constructed and used to investigate uncertainty in the design, production and usage of a mechanical system, see Figure 1a).  [20], c) Sensory rod used for monitoring loads and states In order to gain insights from the system, it is necessary to equip load-bearing structures with sensors and periphery. A novel manufacturing process aims to integrate sensors into the inside of metallic structure, Figure 1b). The sensor is used during its usage phase to measure the axial forces in the rods. By this the conditions of the whole system and the load history can be analyzed, Figure 1c). To qualify the sensor for this measurement task, it has to be designed to show a high sensitivity to axial loads and a low resistance to all other loads. An axial force sensor should not give a signal when applying bending or transverse forces. However, every sensor also has a small sensitivity towards these unintended loads, which is called cross-talk [18]. When designing force sensors, these parasitic effects have to be compensated and minimized as far as possible. For this, numerical models are often used [32]. However, the established design approach does not take into account the inevitable manufacturing uncertainty occurring during production. Therefore, it can be expected that an optimum solution for an ideal sensor element is not the best solution for actual sensors produced under real conditions. For this reason, the task of this investigation is to find an optimum geometry of a sensor element subjected to uncertainty. Uncertain parameters occurring during the production phase are anticipated and used as an input for the robust optimization task.
To find an appropriate shape for the sensor body we formulate a nonlinear PDE-constrained optimization problem where both the objective and the constraints are nonconvex uncertain functions in the sense that they depend on uncertain parameters. Many examples exist in the literature which demonstrate that simply ignoring the uncertain nature of the problem results in solutions that are worthless in practice. In fact, a solution thus found may be highly suboptimal or, worse, infeasible even under slight perturbations of the parameter [3,2,4]. The two most prominent ways to deal with uncertainty appropriately in mathematical programming are stochastic programming [7] and robust optimization [1].
Stochastic programming assumes that the uncertainty is random and a (sufficiently accurate) probabilistic description of it is available. The uncertainty is incorporated into the optimization problem by probabilistic or risk measures -for example, the expected value of some performance measure is optimized. Uncertainty in the constraints leads to probabilistic, or chance, constraints that are appropriate only if violations of the constraints can be tolerated. Moreover, they often are difficult to trea [7, pp. 103 sqq.].
Robust optimization assumes that the uncertain parameters are restricted to a compact uncertainty set with often simple structure, and that the constraints are hard in the sense that violations must not occur [1]. The original problem is reformulated by taking the worst-case behavior into account. The resulting optimization problem is a bilevel or semi-infinite problem (depending on the formulation) called the robust counterpart. The robust counterpart reduces the feasible set to those designs that are robust feasible, i.e., feasible for all possible realizations of the uncertain parameter.
The optimization problem that we will state models uncertain placement of strain gauge sensors. The sensors are applied like an adhesive strip on the surface of a body to measure length variation due to external forces. The signals of the sensors can only be trusted within a certain range specified by the manufacturer, which leads to hard constraints and makes robust optimization the appropriate choice for our problem.
A major part of the literature on robust optimization addresses optimization problems whose robust counterpart can be reformulated as a computationally tractable single-level problem. Examples are linear programs, second-order cone programs, quadratically constrained quadratic programs and semidefinite programs [4]. General nonlinear programs, on the other hand, are much less studied.
An approach known from semi-infinite programming that is applicable to some nonlinear robust optimization problems is the local reduction approach [29,34,35]. Here, the continuous uncertainty set is replaced by the (assumed to be) discrete and finite subset of the (assumed to be) isolated local maximizers of the constraint functions. In general, that subset is found by performing many local optimizations for a fine grid of starting points. The local maximizers so found are tracked by their necessary optimality conditions of first order, which results in a problem of the class of mathematical programs with complementarity constraints (MPCC). A related method is the adversarial approach [19]. The local reduction approach has been applied to an optimal control problem constrained by an ordinary differential equation in [16], see also [23].
Bertsimas et al. compute directions pointing away from all worst-case realizers and show that such directions are descent directions for the (unconstrained) robust counterpart, which is a minimax problem [6]. An extension to constrained optimization has also been investigated [5]. The set of all maximizers is computed in a slightly simpler way than in the local reduction approach.
When both the size of the uncertainty set and the order of nonlinearity of the uncertain functions are moderate. By replacing the uncertain functions with low-order Taylor approximations, the inner maximization problems become computationally tractable and can be solved globally. This method has been investigated with first-order approximations in [36,21,15] and with second-order approximations in [31,26]. The resulting problem is either stated as an MPCC or, if issues due to nondifferentiability are not expected, as a nonlinear program with explicit worst-case functions. Unlike the aforementioned approaches, the second-order MPCC methods cope well with nonisolated worst cases, because they track the global solution by a necessary and sufficient optimality system that is often differentiable.
In this article, we consider a nonsmooth formulation of the robust counterpart that utilizes the value, or worst-case, functions of the lower-level problems. While unfortunately not differentiable everywhere, these functions are nevertheless quite regular and, in fact, differentiable almost everywhere. We will point out in the next section that they are locally Lipschitz continuous with readily available subgradients. The resulting optimization problem is solved with a BFGS-SQP method adapted for nonsmooth problems.

Model Description
Physical description. To optimize the sensor geometry inside the tubular structure, it is required to consider also the adjacent part of the outer tube. This steel tube (E = 210 GPa, d = 27 mm; w = 1.5 mm; length = 130 mm) is fixed at the bottom and loaded with different loads at the upper surface of the free tube end. A tensile force f ten of 12.5 kN is applied in z-direction in the first step. Subsequently, the transverse forces f j tran of 2 kN are applied in the second to the forth step one by one. By this, the cross-talk of the sensor is identified. The steel sensor (E = 210 GPa, l = 42 mm) is arranged in the middle of the tube. The steel tube and the sensor element are depicted in Figure 2. Due to the joining process the end caps are chamfered and the sensor element is pre-tensioned between the undercuts inside the tube. However, in this investigation, this joint is considered as an ideal connection, which can also be achieved in practice if the pre-tension is dimensioned properly [25]. Strain gauge sensors are virtually attached on the outer surface of the sensor element in the middle of the bars. The maximum strains under nominal load (f ten ) are constrained to 4 000 µm/m, as the sensor element is only meant to deform elastically. However, in order to guarantee measurability of the signal, the virtual strain gauge sensors have to experience an elongation of at least 800 µm/m under nominal load (f ten ). Otherwise the signal/noise ratio becomes too low. Each strain gauge is simplified as a line element between two nodes. The relative length change of the line element is corresponding to the signal of each individual strain

134
Uncertainty in Mechanical Engineering III gauge. In practice, the strain gauges are amplified by wiring the variable resistances to a Wheatstone bridge. By this, the signals are amplified and at the same time disturbances are compensated to a certain extent [18]. Mathematically, the strain gauges are summed up, which in theory, compensates all disturbing cross-talk effects caused by bending. However, due to inaccurate bonding of the strain gauge sensors, inhomogeneous material or geometrical deviations, a perfect compensation is never achieved in reality. To consider the uncertainty in the design of the sensor element, the positions of the strain gauge sensors p j i ,i = 1, . . . , 4,j = 1, 2, may vary, details given in the next section. The parametric optimization problem. We proceed with a description of the optimization problem that formalizes the task of finding a properly shaped sensor element. We consider a family of three-dimensional domains { parametrized by a vector µ ∈ X, which is our design variable. For a given design µ ∈ X, the domain Ω(µ) describes the shape of the steel tube and the integrated sensor element. It will be useful to consider a decomposition Ω(µ) = Ω T ∪ Ω S (µ) of the domain into the steel tube Ω T ⊆ R 3 and the sensor element Ω S (µ) ⊆ R 3 , as well as a decomposition Ω S (µ) = Ω SI ∪ Ω SD (µ) of the sensor element into the parts Ω SD (µ) ⊆ R 3 that do and the parts Ω SI ⊆ R 3 that do not depend on the design. 3 Overall, the body decomposes into We denote by Γ(µ) ⊆ R 3 the surface of Ω(µ), and define an analogous decomposition where Γ 0 ⊆ R 3 is the bottom of the tube, which is fixed. Our objective is to find a design such that the sensor element is sensitive to one set of forces and insensitive to a second set of forces, details given below.
The response of an elastic body to small or moderate forces can be described with sufficient accuracy by the equations of linear elasticity, an elliptic linear PDE given by div σ(y) = 0, in Ω(µ), Here, div σ(y) denotes the divergence of the Cauchy stress tensor, n : Γ(µ) → R 3 is the outer unit normal on Γ(µ), and f : Γ T → R 3 is the surface force acting on the body. The unique solution y : Ω(µ) → R 3 defines the displacement of the body in reaction to the surface load. For brevity, we condense the equations of linear elasticity into the single state equation The signal of a virtual strain gauge sensor (VSGS) can be defined in terms of the displacement by where p 1 , p 2 are the two end points of the line element that models the strain gauge sensor, and ∥·∥ denotes the Euclidean norm of R 3 . A negative value indicates that the sensor body to which the VSGS is attached is clinched locally, and a positive value indicates that it is stretched locally. Finally, we define -formally, without derivation -the reduced function s(µ, p 1 , p 2 ; f ) = S µ (y, p 1 , p 2 ), where y solves e µ (y, f ) = 0.
The mapping s is continuously differentiable with respect to µ and continuous with respect to p 1 and p 2 . It can be defined in a rigorous way by a transformation method, see [13, pp.287 sqq.]. We consider four VSGSs defined by four pairs (p i 1 , p i 2 ) ∈ Γ SI × Γ SI , i = 1, . . . , 4, of points, which we pack into a single variable p = (p 1 1 , . . . , p 4 2 ) ∈ P ⊆ (R 3 ) 8 for brevity. We are interested in the signals of the VSGSs in four different load scenarios. The first load scenario is a tensile force f ten : Γ L → R 3 , and the other three scenarios are transverse forces f j tra : Γ L → R 3 , j = 1, 2, 3. We seek a design that makes the sensor element sensitive to the tensile forces but insensitive to the transverse forces. An appropriate choice for the objective function therefore is the sums of the quotients of sensor signals Given positions p ∈ P of the four VSGS, the problem of finding a suitable design is given by where b l , b u ∈ R are scalar values for the lower and upper bounds. The intended positions of the four VSGSs are depicted as black vertical lines between the dark gray rectangles on the crosshatched stripes in Figure 3. We denote that position byp ∈ P , which we also call the nominal parameter. Likewise, we call the optimization problem P(p) the nominal problem. The robust optimization problem. In practice, the VSGSs cannot be mounted exactly to a desired position -in our case, top -but small deviations have to be expected. This means, in problem (P(p)), we cannot expect the parameter p ∈ P to be exactlyp. Instead, we follow the fundamental assumption of robust optimization [1] that the parameter p is restricted to a compact uncertainty set U ⊆ P centered aroundp, but otherwise is unknown. In our case, the uncertainty set is given by where each of the sets U i k ⊆ Γ SI , i = 1, . . . , 4, k = 1, 2, is a rectangle in the angular-vertical plane in cylindrical coordinates centered aroundp i k ∈ Γ SI . If we let (ρ i k ,θ i k ,z i k ) ∈ Γ SI be a representation of the pointp i k in cylindrical coordinates, we set with positive scalars ∆ θ and ∆ z . The uncertainty sets are depicted as magenta rectangles in Figure 3.
In particular, we assume no knowledge about a probability distribution of the parameter. A consequence of this assumption is that we have to make the design decision for µ ∈ X before we actually know the precise value for p, though clearly, the value of p affects optimality and feasibility substantially. An optimal solution, sayμ * ∈ X, of the nominal problem P(p) may be seriously suboptimal or even infeasible for another problem P(p), p ∈ P . This is further discussed, for example, in [3,2].
The robust optimization paradigm deals with this issue with the notion of guaranteed objective and constraint function values. This means that a design µ is considered feasible only if it is feasible for all problems P(p), p ∈ P , then said to be robust feasible. Likewise, the objective function value associated to µ is the worst value of q(µ, p) for any p ∈ P , because this is the only value that can be guaranteed. The resulting optimization problem is the robust counterpart min µ max p∈U q(µ, p) Usually, the robust counterpart has a smaller feasible set than the nominal problem. It may even happen that the robust feasible set is empty, in which case the robustness requirement is too demanding and cannot be satisfied. A shorter but equivalent formulation of the robust counterpart uses the worst-case functions of the objective as well as the constraint functions. These are defined as where the identity min f = − max(−f ) for a continuous function over a compact set has been used. The robust counterpart then can be written as Generally, the solution of an optimization problem is not a differentiable function of the data, and consequently, there may be points where the worst-case functions fail to be differentiable. Nonetheless, the worst-case functions are quite regular, as the following proposition shows.

The (Clarke) generalized gradient is given by
where conv denotes the convex hull.

If M (µ) = {p * } is a singleton, then the generalized gradient ∂Φ
4. Φ is differentiable at almost every µ ∈ X, i.e., the set of points where Φ fails to be differentiable has Lebesque measure zero.
Proof. The first two properties follow from Theorem 2.1 of [9], the third from Proposition 1.13 of that reference, and the last property is due to Rademacher's theorem, see Theorem VIII.3 of [33].
Remark. In many practically relevant applications, the evaluation of a worst-case function -or equivalently, determining an element of M (µ) -is a difficult nonconvex optimization problem on its own that needs to be solved globally. When this is not computationally tractable, approximation techniques have been used [36,21,15,31,26], which, however, require that the uncertain functions are at least once or twice continuously differentiable with respect to the uncertain parameters. In our case, on the other hand, the rather simple form of the uncertain function S µ in (7) enables us to evaluate the functions Φ i and the respective sets of maximizers cheaply, even though S µ are only continuous but not differentiable with respect to the uncertain parameters; the nondifferentiability is due to the nondifferentiability of the solution y across element boundaries in the finite-element discretized setting. Since the number of surface elements in the uncertainty sets is quite small (around 50 each), we can maximize S µ on combinations (around 50 2 = 2 500 per VSGS) of elements. As a consequence, we can compute exact subgradients of Φ i .
Due to Proposition 1, the robust counterpart (RC) is eligible to a variety of optimization methods for nonsmooth problems. Popular among them are bundle methods, see [30] for the case without nonsmooth constraints, [17] for convex nonsmooth objective and constraints, and [28] for nonconvex nonsmooth objective and constraints. Bundle methods require an arbitrary element of the generalized gradient of each nonsmooth function in every iteration. Another approach uses gradient sampling to obtain descent directions, see [8] for the unconstrained and [11] for the constrained case. Sampling methods have appealing theoretic properties even in the nonconvex case [24]. However, they require O(n) gradient evaluations in every iteration, n corresponding to the number of design variables. Lately, quasi-Newton methods originally for smooth problems have been studied for nonsmooth problems, see [27,10]. The authors demonstrate that SQP methods with BFGS updates were more successful and more efficient than the sampling methods on a large test set, despite the theoretical convergence results favoring the sampling methods. Motivated by these promising results, we use the BFGS-SQP method presented in [10] to solve the robust counterpart (RC). The results are presented in the next section.

138
Uncertainty in Mechanical Engineering III

Results
We have solved both the nominal problem P(p) and the robust counterpart (12) and compare the results in this section. We followed the first discretize, then optimize approach, where the PDE given by eqs. (2) to (5) is discretized before the optimization, see Section 3.2.3 of [22]. Since all functions in the nominal problem are differentiable with respect to µ, it can be solved with gradient-based methods. We used the fmincon command of Matlab 2017b, configured to use an ordinary BFGS-SQP method for smooth problems. For the nonsmooth robust counterpart, we used the nonsmooth BFGS-SQP method of [10], a Matlab package of which is publicly available under the name GRANSO 4 ; we used version 1.5.2 of GRANSO. We also applied fmincon to the robust counterpart, to see whether the nonsmoothness actually causes troubles or could simply be ignored. We used the standard termination criteria for both methods. For fmincon, this is a KKT residual of at most 1 · 10 −6 . For GRANSO, this is a solution norm of at most 1 · 10 −6 of a quadratic problem that uses gradient information of several recent and nearby iterations to reduce the penalty functions, see [11] for details. The shape of the vertical borders of the rectangular holes is described with one cubic Bézier curve with four control points each. The first and the fourth control point are two corners (one above the other) of the rectangles and do not change. The middle two control points of each Bézier curve in the angular-vertical plane act as our design variables, which accounts for a total of 32 design variables (eight borders with four degrees of freedom each). We have employed box constraints a ≤ µ i ≤ b for each design variable in order to keep the mesh intact and also to ensure the strips between the rectangles do not become too thin. Since GRANSO does not support box constraints, we used a nonlinear periodic remapping of the variables to the feasible domain for all runs. The designs of the optimal solutions are shown in Figure 4, Figure 6 and Figure 5. Quantitative information is given in Table 1. The respective optimal solution denoted by µ * , we display 1. q(µ * ,p), the objective value in the undisturbed/nominal case, 2. max q(µ * ,p), the objective value in the worst case, 3. ∥c + (µ * ,p)∥, the constraint violation in the undisturbed/nominal case, 4. max∥c + (µ * ,p)∥, the constraint violation in the worst case, 5. the number of iterations, and 6. the number of objective function and constraint evaluations. The constraint values are normalized to 1, so that a violation of 0.7369 is considerable. The starting point for the nominal problem was the situation depicted in Figure 3 with rectangular holes. For the robust counterparts, the starting points were the solution of the nominal problem. The smooth method fmincon did not converge for the robust counterpart. Instead, it terminated because the step size in the line search had become too small. The terminal solution was still not robust feasible, and had a worse worst-case objective value than the nominal solution. On the other hand, GRANSO terminated regularly and found an optimal solution of the robust counterpart that is indeed robust feasible, and improved the worst-case objective value by about 4.3 % at the cost of a nominal objective value which is about 2.7 % worse.
The number of iterations appears to be quite large. We observed, however, that many iterations could be saved by reducing the termination criteria to 1 · 10 −4 without losing anything significant in the quality of the solution. For example, after only 99 iterations, GRANSO had found a solution that is robust feasible with an objective value less than 0.03 % worse than the final solution found after 197

140
Uncertainty in Mechanical Engineering III iterations, but the stationarity measure was just below 1 · 10 −4 . Less drastic, fmincon found a feasible solution with basically the same objective function value as the final solution after 33 (instead of 41) iterations, but the stationarity measure had, too, just fallen below 1 · 10 −4 .

Discussion
While the numbers of Table 1 show the success of the robust optimization in terms of objective and constraint value, we also wish to give an interpretation of the geometry. Two conclusions can be drawn from the numerical investigation: 1. The applied transverse forces cause bending moments within the sensor. This leads to tensile strains at one bar and compression strains at the opposing bar. Therefore, in both optimization problems, the solver tends to solutions in which opposing bars have a similar geometry. This symmetry can be identified to some extent in the nominal solution and explicitly in the robust solution. Even though there was no symmetry constraint given in the model, the robust solution shows a perfect symmetry for all four bars.
2. The robust geometry consists of equally thinned out bars in the region of the VSGS. By reducing the cross-section in this region of the bars, the occurring strains under nominal force f ten increase. As the solver has the task to increase the sensitivity of the sensor under tensile forces, the proposed solutions appears very plausible.
An advantage of our nonsmooth approach to solve the robust counterpart (RC) which we have not mentioned so far is that the constraints can be reduced to a single constraint by taking the maximum of the worst cases and then replacing Φ i (µ) ≤ 0, i = 1, . . . , 8, by the single constraint Ψ(µ) ≤ 0. Due to a result analogous to Proposition 1 the subgradients of Ψ can be obtained via the subgradients of the active elements. For some solvers, this can help reduce the computational effort of each iteration if some constraints are usually inactive, as one can expect to be the case when the constraints consist of both lower and upper bounds. Our approach allows to include different types of uncertainty in the optimization program. For example, we could consider the magnitude and/or angle of the external forces f ten or f tra uncertain. Moreover, manufacturing tolerances in the shape of the holes of the sensor element could be addressed. In both cases, the approximation techniques of first [15] and second [31,26] order would allow to incorporate approximated worst-case functions in our framework.