1 ? Massachusetts Institute of Technology - Prof. de Weck and Prof. Willcox Engineering Systems Division and Dept. of Aeronautics and Astronautics Multidisciplinary System Design Optimization (MSDO) Sensitivity Analysis Lecture 8 1 March 2004 Olivier de Weck Karen Willcox 2 ? Massachusetts Institute of Technology - Prof. de Weck and Prof. Willcox Engineering Systems Division and Dept. of Aeronautics and Astronautics Today’s Topics Sensitivity Analysis – effect of changing design variables – effect of changing parameters – effect of changing constraints Gradient calculation methods – Analytical and Symbolic – Finite difference – Adjoint methods – Automatic differentiation 3 ? Massachusetts Institute of Technology - Prof. de Weck and Prof. Willcox Engineering Systems Division and Dept. of Aeronautics and Astronautics Standard Problem Definition 1 2 min ( ) s.t. ( ) 0 1,.., ( ) 0 1,.., 1,.., j k u iii J gjm hkm xxxi n d dd x x x A For now, we consider a single objective function, J(x). There are n design variables, and a total of m constraints (m=m 1 +m 2 ). The bounds are known as side constraints. 4 ? Massachusetts Institute of Technology - Prof. de Weck and Prof. Willcox Engineering Systems Division and Dept. of Aeronautics and Astronautics Sensitivity Analysis Sensitivity analysis is a key capability aside from the optimization algorithms we discussed. Sensitivity analysis is key to understanding which design variables, constraints, and parameters are important drivers for the optimum solution x*. The process is NOT finished once a solution x* has been found. A sensitivity analysis is part of post- processing. Sensitivity/Gradient information is also needed by: – gradient search algorithms – isoperformance/goal programming – robust design 5 ? Massachusetts Institute of Technology - Prof. de Weck and Prof. Willcox Engineering Systems Division and Dept. of Aeronautics and Astronautics Sensitivity Analysis How sensitive is the “optimal” solution J* to changes or perturbations of the design variables x*? How sensitive is the “optimal” solution x* to changes in the constraints g(x), h(x) and fixed parameters p ? 6 ? Massachusetts Institute of Technology - Prof. de Weck and Prof. Willcox Engineering Systems Division and Dept. of Aeronautics and Astronautics Sensitivity Analysis: Aircraft Questions for aircraft design: How does my solution change if I change the cruise altitude? change the cruise speed? change the range? change material properties? relax the constraint on payload? .. 7 ? Massachusetts Institute of Technology - Prof. de Weck and Prof. Willcox Engineering Systems Division and Dept. of Aeronautics and Astronautics Sensitivity Analysis Questions for spacecraft design: How does my solution change if I change the orbital altitude? change the transmission frequency? change the specific impulse of the propellant? change launch vehicle? Change desired mission lifetime? .. 8 ? Massachusetts Institute of Technology - Prof. de Weck and Prof. Willcox Engineering Systems Division and Dept. of Aeronautics and Astronautics Gradient Vector – single objective “How does the objective function J value change as we change elements of the design vector x?” 1 2 n J x J x J x w ao ?? w ?? w?? ?? w ? ?? ?? ?? w ?? ?? w ?? J # Compute partial derivatives of J with respect to x i i J x w w ?J Gradient vector points normal to the tangent hyperplane of J(x) 1 x 2 x 3 x 9 ? Massachusetts Institute of Technology - Prof. de Weck and Prof. Willcox Engineering Systems Division and Dept. of Aeronautics and Astronautics Geometry of Gradient vector (2D) 0 0.5 1 1.5 2 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 x 1 x 2 Contour plot 3 .1 3 . 1 3 .1 3 . 2 5 3 . 2 5 3 . 2 5 3. 25 3 . 2 5 3 .5 3 .5 3.5 3 . 5 4 4 4 4 5 5 Example function: 12 1 2 12 1 ,Jxx x x xx  ? 2 11 2 21 1 1 1 1 J xxx J J xxx w aoa o  ??? ? w ??? ? ? w??? ?  ??? ? w ??? ? Gradient normal to contours 10 ? Massachusetts Institute of Technology - Prof. de Weck and Prof. Willcox Engineering Systems Division and Dept. of Aeronautics and Astronautics Geometry of Gradient vector (3D) 222 123 Jxxx  1 2 3 2 2 2 x Jx x ao ?? ? ?? ?? ?? increasing values of J 1 x 2 x 3 x Tangent plane 123 22260xxx > @ 111 T o x > @ 222 o T J? x J=3 Example Gradient vector points to larger values of J 11 ? Massachusetts Institute of Technology - Prof. de Weck and Prof. Willcox Engineering Systems Division and Dept. of Aeronautics and Astronautics Taylor Series Expansion , kk n Jwherxx6\6Taylor Series Expansion of Objective Function Tangential hyperplane at x o Effect of curvature (2nd derivative) at x o 000 0 0 1 () () ()( ) ( )()( )H.O.T. 2 T T JJ J ao ?   ?? 0 xx xxx xxHxxx first order term second order term 12 ? Massachusetts Institute of Technology - Prof. de Weck and Prof. Willcox Engineering Systems Division and Dept. of Aeronautics and Astronautics Jacobian Matrix – multiple objectives If there is more than one objective function, i.e. if we have a gradient vector for each J i , arrange them columnwise and get Jacobian matrix: 12 11 1 12 22 2 12 z z z nn n J JJ xx x J JJ xx x J JJ xx x ww w ao ?? ww w ww w ?? ww w ? ?? ww w ww w ?? J " " ##%# " 1 2 z J J J ao ?? ?? ?? ?? ?? J # n x z z x 1 13 ? Massachusetts Institute of Technology - Prof. de Weck and Prof. Willcox Engineering Systems Division and Dept. of Aeronautics and Astronautics Normalization In order to compare sensitivities from different design variables in terms of their relative sensitivity it is necessary to normalize: i J x w w o x “raw” - unnormalized sensitivity = partial derivative evaluated at point x i,o , () io ii i x JJ J xx J x 'w ? o o x x  Normalized sensitivity captures relative sensitivity ~ % change in objective per % change in design variable Important for comparing effect between design variables 14 ? Massachusetts Institute of Technology - Prof. de Weck and Prof. Willcox Engineering Systems Division and Dept. of Aeronautics and Astronautics Example: Dairy Farm Problem With respect to which design variable is the objective most sensitive? “Dairy Farm” sample problem L R N L – Length = 100 [m] N - # of cows = 10 R – Radius = 50 [m] fence 2 2 22 100 / ALRR FL R M AN S S   ? C f FnN INMm PIC ?? ??  Parameters: f=100$/m n=2000$/cow m=2$/liter x o Assume that we are not at the optimal point x* ! COW COW COW 15 ? Massachusetts Institute of Technology - Prof. de Weck and Prof. Willcox Engineering Systems Division and Dept. of Aeronautics and Astronautics Dairy Farm Sensitivity Compute objective at x o Then compute raw sensitivities Normalize Show graphically (optional) 36.6 2225.4 588.4 P L P J N P R w ao ?? w ao ?? w ?? ?? ? ?? w ?? ?? ?? w ?? ?? w ?? 100 36.6 13092 0.28 10 2225.4 1.7 ( ) 13092 2.25 50 588.4 13092 o o JJ J ao ? ?? ao ?? ?? ? ? ?? ?? ?? ?? ?? x x ( ) 13092 o J x Dairy Farm Normalized Sensitivities 0 0.5 1 1.5 2 2.5 L N R Design Variable 16 ? Massachusetts Institute of Technology - Prof. de Weck and Prof. Willcox Engineering Systems Division and Dept. of Aeronautics and Astronautics Realistic Example: Spacecraft X Y Z What are the design variables that are “drivers” of system performance ? 012 meters Spacecraft CAD model -60 -40 -20 0 20 40 60 -60 -40 -20 0 20 40 60 Centroid X [Pm] Centroi d Y [ P m] J 2 = Centroid Jitter on Focal Plane [RSS LOS] T=5 sec 14.97 Pm 1 pixel Requirement: J 2,req =5 Pm NASA Nexus Spacecraft Concept Finite Element Model Simulation “x”-domain “J”-domain 17 ? Massachusetts Institute of Technology - Prof. de Weck and Prof. Willcox Engineering Systems Division and Dept. of Aeronautics and Astronautics Graphical Representation Graphical Representation of Jacobian evaluated at design x o , normalized for comparison. -0.5 0 0.5 1 1.5 Kcf Kc fca Mgs QE Ro lambda zeta I_propt I_ss t_sp K_zpet m_bus K_rISO K_yPM m_SM Tgs Sst Srg Tst Qc fc Ud Us Ru J1: Norm Sensitivities: RMMS WFE Desi gn Variables o1,o w 1 wx/J *J/x analytical finite difference -0.5 0 0.5 1 1.5 Kcf Kc fca Mgs QE Ro lambda zeta I_propt I_ss t_sp K_zpet m_bus K_rISO K_yPM m_SM Tgs Sst Srg Tst Qc fc Ud Us Ru J2: Norm Sensitivities: RSS LOS wwx o /J 2,o *J 2 / x disturbance variables structural variables control variables optics variables 12 0 12 uu o cf cf J J RR J J J J KK aoww ?? ww ?? ? ww ?? ww ?? x "" J1: RMMS WFE most sensitive to: Ru - upper wheel speed limit [RPM] Sst - star tracker noise 1V [asec] K_rISO - isolator joint stiffness [Nm/rad] K_zpet - deploy petal stiffness [N/m] J2: RSS LOS most sensitive to: Ud - dynamic wheel imbalance [gcm 2 ] K_rISO - isolator joint stiffness [Nm/rad] zeta - proportional damping ratio [-] Mgs - guide star magnitude [mag] Kcf - FSM controller gain [-] 18 ? Massachusetts Institute of Technology - Prof. de Weck and Prof. Willcox Engineering Systems Division and Dept. of Aeronautics and Astronautics Analytical Sensitivities If the objective function is known in closed form, we can often compute the gradient vector(s) in closed form (analytically, symbolically): Example: 12 1 2 12 1 ,Jxx x x xx  ? 2 11 2 21 1 1 1 1 J xxx J J xxx w aoa o  ??? ? w ??? ? ? w??? ?  ??? ? w ??? ? Example x 1 = x 2 =1 J(1,1)=3 0 (1,1) 0 J ao ? ?? ?? Minimum Analytical Gradient: For complex systems analytical gradients are rarely available 19 ? Massachusetts Institute of Technology - Prof. de Weck and Prof. Willcox Engineering Systems Division and Dept. of Aeronautics and Astronautics Symbolic Differentiation Use symbolic mathematics programs E.g. Matlab,Maple, Mathematica EDU? syms x1 x2 EDU? J=x1+x2+1/(x1*x2); EDU? dJdx1=diff(J,x1) dJdx1 =1-1/x1^2/x2 EDU? dJdx2=diff(J,x2) dJdx2 = 1-1/x1/x2^2 construct a symbolic object difference operator 20 ? Massachusetts Institute of Technology - Prof. de Weck and Prof. Willcox Engineering Systems Division and Dept. of Aeronautics and Astronautics Finite Differences (I) Taylor Series expansion Neglect second order and H.O.T. Solve for gradient vector 2 '''2 2 oooo x fx x fx xf x f x O x ' ' '   ' x x o x+'x x-'x 'x 'x f Function of a single variable f(x) ' oo o fx x fx fx Ox x '  ' ' Forward Difference Approximation to the derivative '' 2 oo x Ox f x xx ] ] ' ' dd' Truncation Error 0,xx'! '?21 ? Massachusetts Institute of Technology - Prof. de Weck and Prof. Willcox Engineering Systems Division and Dept. of Aeronautics and Astronautics Finite Differences (II) 1 11 111 1 111 1 1 oo o o Jx Jx Jx x Jx J J xxx x x ' w' | w ' ' J(x) 1 1 x 1 o x 1 1 Jx 1 o Jx true, analytical sensitivity finite difference approximation 1 x' J' 1 1 x 1 o x - x 1 22 ? Massachusetts Institute of Technology - Prof. de Weck and Prof. Willcox Engineering Systems Division and Dept. of Aeronautics and Astronautics Finite Differencing (III) Take Taylor expansion backwards at o xx' 2 '''2 2 oooo x fx x fx xf x f x O x ' ' '   ' 2 '''2 2 oooo x fx x fx xf x f x O x ' ' '   ' (1) (2) (1)-(2) and solve again for derivative '2 2 oo o fx x fx x fx Ox x '  ' ' ' 2 2'' 6 oo x Ox f xxx ] ] ' ' dd' Truncation Error Central Difference Approximation to the derivative 23 ? Massachusetts Institute of Technology - Prof. de Weck and Prof. Willcox Engineering Systems Division and Dept. of Aeronautics and Astronautics Finite Difference Overview x' Forward Difference '( ) oo o fx x fx fx x '  | ' 1 st derivative 2 nd derivative 2 22 ''( ) ooo o fx x fx x fx fx x '  '  | ' x'x' Central Difference 2 nd derivative 1 st derivative '( ) 2 oo o fx x fx x fx x '  ' | ' 2 2 ''( ) ooo o fx x fx fx x fx x '   ' | ' 24 ? Massachusetts Institute of Technology - Prof. de Weck and Prof. Willcox Engineering Systems Division and Dept. of Aeronautics and Astronautics Errors of Finite Differencing Caution: - Finite Differencing always has errors - very dependent on perturbation size 12 1 2 12 1 ,Jxx x x xx  ? x 1 = x 2 =1 J(1,1)=3 0 (1,1) 0 J ao ? ?? ?? 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 0 1 2 3 4 5 6 7 8 9 10 x 1 'x 1 10 -9 10 -8 10 -7 10 -6 10 -8 10 -7 10 -6 10 -5 Truncation Errors ~ 'x Rounding Errors ~ 1/'x Perturbation Step Size 'x 1 Gradient Error Choice of 'x is critical 25 ? Massachusetts Institute of Technology - Prof. de Weck and Prof. Willcox Engineering Systems Division and Dept. of Aeronautics and Astronautics Perturbation Size 'x Choice Error Analysis Significant digits (Barton 1992) Machine Precision Trial and Error – typical value ~ 0.1-1% 1/2 A xfH'# - Forward difference 1/3 A x fH'# - Central difference (Gill et al. 1981) o x theoretical function computed values ~ x' A H 10 q kk xx  '#? Step size at k-th iteration q-# of digits of machine Precision for real numbers 26 ? Massachusetts Institute of Technology - Prof. de Weck and Prof. Willcox Engineering Systems Division and Dept. of Aeronautics and Astronautics Computational Expense of FD Cost of a single objective function evaluation of J i i F J Cost of gradient vector finite difference approximation for J i for a design vector of length n i nFJ? Cost of Jacobian finite difference approximation with z objective functions i z nFJ?? Example: 6 objectives 30 design variables 1 sec per simcode evaluation 3 min of CPU time for a single Jacobian estimate - expensive ! 27 ? Massachusetts Institute of Technology - Prof. de Weck and Prof. Willcox Engineering Systems Division and Dept. of Aeronautics and Astronautics Automatic Differentiation Mathematical formulas are built from a finite set of basic functions, e.g. sin x, cos x, exp x Take analysis code in C or Fortran Using chain rule, add statements that generate derivatives of the basic functions Tracks numerical values of derivatives, does not track symbolically as discussed before Outputs modified program = original + derivative capability 28 ? Massachusetts Institute of Technology - Prof. de Weck and Prof. Willcox Engineering Systems Division and Dept. of Aeronautics and Astronautics Chain Rule example () () uqs spt quantity of interest First compute Want to take derivatives w.r.t “t” () ds d p t dt dt Store this value numerically Then apply chain rule (()) () du d d ds qst qs dt dt ds dt ? substitute desired sensitivity 29 ? Massachusetts Institute of Technology - Prof. de Weck and Prof. Willcox Engineering Systems Division and Dept. of Aeronautics and Astronautics Chain Rule Coding Example hr = gm*eps/rho-0.5*g1*(ux*ux + vy*vy); h_u[0] = -gm*eps/(rho*rho); h_u[1] = -g1*ux; h_u[2] = -g1*vy; h_u[3] = gm/rho; hi = (di*hr+hl)*d1; hi_u[0] = (di*hr+hl)*d1_u[0] + d1*(di_u[0]*hr+di*h_u[0]); compute hr differentiate: wrt rho wrt ux wrt vy wrt eps 30 ? Massachusetts Institute of Technology - Prof. de Weck and Prof. Willcox Engineering Systems Division and Dept. of Aeronautics and Astronautics Adjoint Methods A way to get gradient information in a computationally efficient way Based on theory from controls Applied extensively in aerodynamic design and optimization For example, in aerodynamic shape design, need objective gradient with respect to shape parameters and with respect to flow parameters – Would be expensive if finite differences are used! Adjoint methods have allowed optimization to be used for complicated, high-fidelity fluids problems. 31 ? Massachusetts Institute of Technology - Prof. de Weck and Prof. Willcox Engineering Systems Division and Dept. of Aeronautics and Astronautics Adjoint Methods Consider where J is the cost function, w contains the N flow variables, and F contains the n shape design variables. (,)JJ wF At an optimum, the variation of the cost function is zero: 0 TT JJ JGGG ww ao ao  ?? ?? ?? ?? wF wF Nu1 nu1 N>>n 32 ? Massachusetts Institute of Technology - Prof. de Weck and Prof. Willcox Engineering Systems Division and Dept. of Aeronautics and Astronautics Adjoint Methods Fluid governing equations: (,) 0R wF 0 RR RGGG ww ao ao  ?? ?? ?? ?? wF wF We can append these constraints to the cost function using a Lagrange multiplier approach: T TT TT JJ RR J JR JR G GGMGG M GMG §·ww ww ao ao ao ao   ¨? ?? ?? ?? ?? ?? ?? ?? ?? ?1 §·§· ww ww ao ao ao ao   ¨?¨? ?? ?? ?? ?? ?? ?? ?? ?? ?1?1 wF wF wF wF wF ww FF 33 ? Massachusetts Institute of Technology - Prof. de Weck and Prof. Willcox Engineering Systems Division and Dept. of Aeronautics and Astronautics Adjoint Methods TT JR JR JGMGMG §·§· ww ww ao ao ao ao   ¨?¨? ?? ?? ?? ?? ?? ?? ?? ?? ?1?1 wF ww FF Choose M to satisfy the adjoint equation: T RJ M ww ao ao ?? ?? ?? ?? ww equivalent to one flow solve T T JR JGMG §· ww ao ao  ¨? ?? ?? ?? ?? ?1 F FF Then does not depend on the number of flow variables total gradient of J 34 ? Massachusetts Institute of Technology - Prof. de Weck and Prof. Willcox Engineering Systems Division and Dept. of Aeronautics and Astronautics Sensitivity Analysis “How does the optimal solution change as we change the problem parameters?” effect on design variables effect on objective function effect on constraints Want to answer this question without having to solve the optimization problem again. Two approaches: – use Kuhn-Tucker conditions – use feasible directions 35 ? Massachusetts Institute of Technology - Prof. de Weck and Prof. Willcox Engineering Systems Division and Dept. of Aeronautics and Astronautics Parameters Parameters p are the fixed assumptions. How sensitive is the optimal solution x* with respect to fixed parameters ? Optimal solution: x* =[ R=106.1m, L=0m, N=17 cows] T Example: “Dairy Farm” sample problem L R N fence Fixed parameters: Parameters: f=100$/m - Cost of fence n=2000$/cow - Cost of a single cow m=2$/liter - Market price of milk How does x* change as parameters change? Maximize Profit COW COW COW 36 ? Massachusetts Institute of Technology - Prof. de Weck and Prof. Willcox Engineering Systems Division and Dept. of Aeronautics and Astronautics Sensitivity Analysis Recall the Kuhn-Tucker conditions. Let us assume that we have M active constraints, which are contained in the vector * ? (*) ( ) 0 ? (*) 0, 0, jj jM j j Jg gjM jM O O ? ?? ? !? | xx x ? ()gx For a small change in a parameter, p, we require that the Kuhn-Tucker conditions remain valid: (KT conditions) 0 d dp 37 ? Massachusetts Institute of Technology - Prof. de Weck and Prof. Willcox Engineering Systems Division and Dept. of Aeronautics and Astronautics Sensitivity Analysis First, let us write out the components of the first equation: ** ? () ()0 jj jM JgO ? ?? | xx ** ? () ()0, 1,.., j j jM ii g J in xx O ? w w  ww | Now differentiate with respect to the parameter p using the chain rule: 1 n i k i dY Y Y x dp p x p www  www | 38 ? Massachusetts Institute of Technology - P Prof. de Weck and Prof. Willcox Engineering Systems Division and Dept. of Aeronautics and Astronautics Sensitivity Analysis ** ? () ()0 j j jM ii g J xx O ? w w  ww | ? (*) 0 j g x differentiate wrt p: 2 2 1 2 2 ? ?? 0 n j k j kjM ik ik jjj jM jM iii g Jx xx xx p gg J xp xp x p O O ? ?? ao w ww  ?? ww ww w ?? www w  ww ww w w || || 1 ?? 0 n jj k k k gg x pxp ww w  www | unknowns are and j i x pp Ow w ww 39 ? Massachusetts Institute of Technology - Prof. de Weck and Prof. Willcox Engineering Systems Division and Dept. of Aeronautics and Astronautics Sensitivity Analysis In matrix form we can write: 0 0 T AB c Bd G G -?-?ao  ???? ?? ˉ?ˉ??? x  n M Mn 2 2 ? j ik j jM ik ik g J A xx xx O ? w w  ww ww | ? j ij i g B x w w 2 2 ? j i jM ii g J c xp xp ? w w  ww ww | ? j j g d p w w 1 2 n x p x p x p G w -? °° w °° w °° °° w ?? °° °° w °° °° w ˉ? x # 1 2 M p p p O O G O w -? °° w °° w °° °° w ?? °° °° w °° °° w ˉ?  # 40 ? Massachusetts Institute of Technology - Prof. de Weck and Prof. Willcox Engineering Systems Division and Dept. of Aeronautics and Astronautics Sensitivity Analysis We solve the system to find Gx and GOO, then the sensitivity of the objective function with respect to p can be found: T dJ J J dp p G w ? w x dJ Jp dp '| ' (first-order approximation) pG'| 'xx To assess the effect of changing a different parameter, we only need to calculate a new RHS in the matrix system. 41 ? Massachusetts Institute of Technology - Prof. de Weck and Prof. Willcox Engineering Systems Division and Dept. of Aeronautics and Astronautics Sensitivity Analysis - Constraints We also need to assess when an active constraint will become inactive and vice versa An active constraint will become inactive when its Lagrange multiplier goes to zero: j jj pp p O OGO w ' ' ' w Find the 'p that makes O j zero: 0 jj pOGO' j j pjM O GO  ' ? This is the amount by which we can change p before the j th constraint becomes inactive (to a first order approximation) 42 ? Massachusetts Institute of Technology - Prof. de Weck and Prof. Willcox Engineering Systems Division and Dept. of Aeronautics and Astronautics Sensitivity Analysis - Constraints An inactive constraint will become active when g j (x) goes to zero: () (*) (*) 0 T jj j gg pgG ao '? ?? xx xx Find the 'p that makes g j zero: (*) (*) j T j g p g G  ' ? x xx for all j not active at x* This is the amount by which we can change p before the j th constraint becomes active (to a first order approximation) If we want to change p by a larger amount, then the problem must be solved again including the new constraint Only valid close to the optimum 43 ? Massachusetts Institute of Technology - Prof. de Weck and Prof. Willcox Engineering Systems Division and Dept. of Aeronautics and Astronautics Lecture Summary Sensitivity analysis – Yields important information about the design space, both as the optimization is proceeding and once the “optimal” solution has been reached. Gradient calculation approaches – Analytical and Symbolic – Finite difference – Automatic Differentiation – Adjoint methods Reading Papalambros – Section 8.2 Computing Derivatives