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) Numerical Optimization II Lecture 7 25 February 2004 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 ? Sequential Linear Programming ? Penalty and Barrier Methods ? Sequential Quadratic Programming ? Mixed Integer Programming 3 ? Massachusetts Institute of Technology - Prof. de Weck and Prof. Willcox Engineering Systems Division and Dept. of Aeronautics and Astronautics Technique Overview Steepest Descent Conjugate Gradient Quasi-Newton Newton Simplex – linear SLP – often not effective SQP – nonlinear, expensive, common in engineering applications Exterior Penalty – nonlinear, discontinuous design spaces Interior Penalty – nonlinear Generalized Reduced Gradient – nonlinear Method of Feasible Directions – nonlinear Mixed Integer Programming UNCONSTRAINED CONSTRAINED 4 ? 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 ≤= == ≤≤ = x x x null 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 ). For now we assume all x i are real and continuous. 5 ? Massachusetts Institute of Technology - Prof. de Weck and Prof. Willcox Engineering Systems Division and Dept. of Aeronautics and Astronautics Optimization Process x 0 , q=0 Calculate ?J(x q ) Calculate S q Perform 1-D search x q = x q-1 + α S q Converged? q=q+1 no yes Done 6 ? Massachusetts Institute of Technology - Prof. de Weck and Prof. Willcox Engineering Systems Division and Dept. of Aeronautics and Astronautics Constrained Optimization Definitions: ? Feasible design: a design that satisfies all constraints ? Infeasible design: a design that violates one or more constraints ? Optimum design: the choice of design variables that minimizes the objective function while satisfying all constraints In general, constrained optimization algorithms try to cast the problem as an unconstrained optimization and then use one of the techniques we looked at in Lecture 6. 7 ? Massachusetts Institute of Technology - Prof. de Weck and Prof. Willcox Engineering Systems Division and Dept. of Aeronautics and Astronautics Linear Programming Most engineering problems of interest are nonlinear ? Can often simplify nonlinear problem by linearization ? LP is often the basis for developing more complicated NLP algorithms Standard LP problem: All LP problems can be converted to this form. = = = == ≥= ∑ ∑ 1 1 min ( ) 1,.., 0 1,..., n ii i n ij i j i i Jcx ax b j m xi n x min ( ) 01,., T i J xi n = = ≥= xcx Ax b 8 ? Massachusetts Institute of Technology - Prof. de Weck and Prof. Willcox Engineering Systems Division and Dept. of Aeronautics and Astronautics Linear Programming To convert inequality constraints to equality constraints, use additional design variables: where x n+1 ≥ 0 x n+1 is called a slack variable e.g. the constraint x 1 +x 2 ≤1 can be written x 1 +x 2 -x 3 =1 x i ≥0, i=1,2,3 1 11 becomes nn jii j jii n j ax b ax x b + == ≤+= ∑∑ if x 3 =0, this constraint is active 9 ? Massachusetts Institute of Technology - Prof. de Weck and Prof. Willcox Engineering Systems Division and Dept. of Aeronautics and Astronautics Simplex Method Solutions at the “vertices” of the design space are called basic feasible solutions. The Simplex algorithm moves from BFS to BFS so that the objective always improves. feasible region basic feasible solution 10 ? Massachusetts Institute of Technology - Prof. de Weck and Prof. Willcox Engineering Systems Division and Dept. of Aeronautics and Astronautics Sequential Linear Programming Consider a general nonlinear problem linearized via first order Taylor series: 00 00 00 min ( ) ( ) ( ) s.t. ( ) ( ) ( ) 0 ( ) ( ) ( ) 0 T T jj j T kk k u ii ii JJ J gg g hh h xx xx δ δ δ δ ≈+? ≈ +? ≤ ≈ +? = ≤+ ≤ xx xx xx xx xx xx null 0 δ = ?xxx where This is an LP problem with the design variables contained in δx. The functions and gradients evaluated at x 0 are constant coefficients. 11 ? Massachusetts Institute of Technology - Prof. de Weck and Prof. Willcox Engineering Systems Division and Dept. of Aeronautics and Astronautics Sequential Linear Programming 1. Initial guess x 0 2. Linearize about x 0 using first-order Taylor series 3. Solve resulting LP to find δx 4. Update: x 1 = x 0 + δx 5. Linearize about x 1 and repeat: x q = x q-1 + δx where δx is the solution of an LP (model linearized about x q-1 ). 12 ? Massachusetts Institute of Technology - Prof. de Weck and Prof. Willcox Engineering Systems Division and Dept. of Aeronautics and Astronautics Sequential Linear Programming ? Linearization approximation is only valid close to x 0 ? Need to restrict size of update δx ? Not considered to be a good method 13 ? Massachusetts Institute of Technology - Prof. de Weck and Prof. Willcox Engineering Systems Division and Dept. of Aeronautics and Astronautics Nonlinear vs. Linear Constraints Linear constraints: ? start from initial feasible point, then all subsequent iterates are feasible ? can construct search direction & step length so that constraints are satisfied ? improvement from x k to x k+1 based on J(x) Nonlinear constraints: ? not straightforward to generate a sequence of feasible iterates ? if feasibility is not maintained, then need to decide whether x k+1 is “better” than x k feasible region 14 ? Massachusetts Institute of Technology - Prof. de Weck and Prof. Willcox Engineering Systems Division and Dept. of Aeronautics and Astronautics Merit Function Is objective reduced? merit function Are constraints satisfied? merit function = f(J(x), g(x), h(x)) examples: penalty function methods, barrier function methods 15 ? Massachusetts Institute of Technology - Prof. de Weck and Prof. Willcox Engineering Systems Division and Dept. of Aeronautics and Astronautics Subproblems ? Many optimization algorithms get to the optimum by generating and solving a sequence of subproblems that are somehow related to the original problem. ? Typically there are two tasks at each iteration: 1. Calculate search direction 2. Calculate the step length ? Sometimes, the initial formulation of a subproblem may be defective i.e. the subproblem has no solution or the solution is unbounded ? A valid subproblem is one for which a solution exists and is well defined ? The option to abandon should be available if the subproblem appears defective ? It is possible that a defective subproblem is an accurate reflection of the original problem having no valid solution 16 ? Massachusetts Institute of Technology - Prof. de Weck and Prof. Willcox Engineering Systems Division and Dept. of Aeronautics and Astronautics Penalty and Barrier Methods General Approach: – minimize objective as unconstrained function – provide penalty to limit constraint violations – magnitude of penalty varies throughout optimization – called sequential unconstrained minimization techniques (SUMT) – create pseudo-objective: (, ) () () pp rJ rPΦ =+xxx J(x) = original objective function P(x) = imposed penalty function r p = scalar multiplier to determine penalty magnitude p = unconstrained minimization number 17 ? Massachusetts Institute of Technology - Prof. de Weck and Prof. Willcox Engineering Systems Division and Dept. of Aeronautics and Astronautics Penalty Methods Four basic methods: (i) Exterior Penalty Function Method (ii) Interior Penalty Function Method (Barrier Function Method) (iii) Log Penalty Function Method (iv) Extended Interior Penalty Function Method Effect of penalty function is to create a local minimum of unconstrained problem “near” x*. Penalty function adds a penalty for infeasibility, barrier function adds a term that prevents iterates from becoming infeasible. 18 ? Massachusetts Institute of Technology - Prof. de Weck and Prof. Willcox Engineering Systems Division and Dept. of Aeronautics and Astronautics Exterior Penalty Function Method (, ) () () pp JPρ ρΦ =+xxx {} [] 12 2 2 11 () max0, () () mm jk jk Pgh == ??=+ ?? ∑∑ xxx ? if all constraints are satisfied, then P(x)=0 ? ρ p = penalty parameter; starts as a small number and increases ?if ρ p is small, Φ(x,ρ p ) is easy to minimize but yields large constraint violations ?if ρ p is large, constraints are all nearly satisfied but optimization problem is numerically ill-conditioned ? if optimization stops before convergence is reached, the design will be infeasible 19 ? Massachusetts Institute of Technology - Prof. de Weck and Prof. Willcox Engineering Systems Division and Dept. of Aeronautics and Astronautics Quadratic Penalty Function [] 12 ? 2 2 11 ? (, ) () () () mm Qp p j k jk Jghρρ == ?? ?? Φ=+ + ?? ?? ?? ∑∑ xx x x ? contains those inequality constraints that are violated at x. ? It can be shown that ? Φ Q (x,ρ p ) is well-defined everywhere ? () j g x lim ( ) p p ρ ρ →∞ =x* x* Algorithm: -choose ρ 0 , set k=0 -find min Φ Q (x,ρ k ) x k * -if not converged, set ρ k+1 > ρ k , k←k+1 and repeat ? Note: x k * could be infeasible wrt the original problem 20 ? Massachusetts Institute of Technology - Prof. de Weck and Prof. Willcox Engineering Systems Division and Dept. of Aeronautics and Astronautics Quadratic Penalty Function Example Φ Q (x, r p ) 10 J(x) 2 min ( ) s.t. 1 0 Jx x x = ?= ρ = 1 * ρ = 10 * ρ = 5 * ρ = 2 * Φ Q (x,ρ)=x 2 + ρ(1-x) 2 x 21 ? Massachusetts Institute of Technology - Prof. de Weck and Prof. Willcox Engineering Systems Division and Dept. of Aeronautics and Astronautics Absolute Value Penalty Function 12 ? 11 ? (, ) () () () mm Ap p j k jk Jghρρ == ?? Φ=+ + ?? ?? ∑∑ xx x x ? () j g x where contains those inequality constraints which are violated at x. ? Φ A has discontinuous derivatives at points where or are zero. ? The crucial difference between Φ Q and Φ A is that we do not require ρ p →∞ for x* to be a minimum of Φ A , so we can avoid ill-conditioning ? Instead, for some threshold , x* is a minimum of Φ A for any ? Sometimes called exact penalty functions ? () j g x () k h x p ρ p p ρ ρ> 22 ? Massachusetts Institute of Technology - Prof. de Weck and Prof. Willcox Engineering Systems Division and Dept. of Aeronautics and Astronautics Absolute Value Penalty Function Example Φ A (x,ρ) 1 0 J(x) ρ = 2 ρ = 5 ρ = 10 x 2 min ( ) s.t. 1 0 Jx x x = ?= * ρ = 1 ρ = 1.5 * Φ Α (x,ρ)=x 2 + ρ |1-x| 23 ? Massachusetts Institute of Technology - Prof. de Weck and Prof. Willcox Engineering Systems Division and Dept. of Aeronautics and Astronautics Interior Penalty Function Method (Barrier Function Method) 1 ? 1 1 () ? () m j j P g = ? = ∑ x x P j (x) g j (x) [] 12 ? 2 11 1 (, , ) () () ? () mm pp p p k jk j rJr h g ρρ == ? Φ=+ + ∑∑ xx x x ? r p = barrier parameter; starts as a large positive number and decreases ? barrier function for inequality constraints only ? sequence of improving feasible designs ? Φ(x,r p ,ρ p ) discontinuous at constraint boundaries 24 ? Massachusetts Institute of Technology - Prof. de Weck and Prof. Willcox Engineering Systems Division and Dept. of Aeronautics and Astronautics Barrier Function Method 1 1 () ln () m j j Pg = ? ?=?? ? ? ∑ xx ? often better numerically conditioned than ? penalty function has a positive singularity at the boundary of the feasible region ? penalty function is undefined for g i (x)>0 ? 1 () j g ? x () 1 ? 1 ? (, , ) () ln () m Lpp p j j rJr gρ = Φ=?? ∑ xx x 0 p r lim ( ) p r → =x* x* 25 ? Massachusetts Institute of Technology - Prof. de Weck and Prof. Willcox Engineering Systems Division and Dept. of Aeronautics and Astronautics Barrier Function Example J(x) r = 5 r = 1 r = 0.5 r = 0.1 Φ L (x,r) 2 min ( ) s.t. 1 0 Jx x x = ?≥ 1 0 10 Φ L (x,r)=x 2 -r ln(x-1) x 26 ? Massachusetts Institute of Technology - Prof. de Weck and Prof. Willcox Engineering Systems Division and Dept. of Aeronautics and Astronautics Extended Interior Penalty Function Method Combine features of interior/exterior methods 1 1 2 () () 1 where ( ) if ( ) () 2() if ( ) m j j jj j j j Pg gg g g g ε ε ε ε = = ? =≤ ? = ?> ∑ xx xx x x x null null -linear extended penalty function -ε is a small negative number, and marks transition from interior penalty to extended penalty -ε must be chosen so that Φ has positive slope at constraint boundary -can also use quadratic extended and variable penalty functions 27 ? Massachusetts Institute of Technology - Prof. de Weck and Prof. Willcox Engineering Systems Division and Dept. of Aeronautics and Astronautics Sequential Quadratic Programming ? Create a quadratic approximation to the Lagrangian ? Create linear approximations to the constraints ? Solve the quadratic problem to find the search direction, S ? Perform the 1-D search ? Update the approximation to the Lagrangian 28 ? Massachusetts Institute of Technology - Prof. de Weck and Prof. Willcox Engineering Systems Division and Dept. of Aeronautics and Astronautics Sequential Quadratic Programming Create a subproblem with quadratic objective function: T 1 min ( ) ( ) ( ) 2 kk kkT QJ J=+? +Sx xSSBS and linear constraints T 1 T 2 s.t. () ()0 1,, () ()0 1,, kk k jj kk k ggjm hhk ?+≤= ?+== xS x xS x … … ? Design variables are the components of S ? B=I initially, then is updated to approximate H (as in quasi-Newton) 29 ? Massachusetts Institute of Technology - Prof. de Weck and Prof. Willcox Engineering Systems Division and Dept. of Aeronautics and Astronautics Incompatible Constraints ? The constraints of the subproblem can be incompatible, even if the original problem has a well-posed solution ? For example, two linearized constraints could be linearly dependent ? This is a common occurrence in practice ? Likelihood of incompatible constraints reduced by allowing flexibility in RHS (e.g. allow scaling factors in front of g j (x k ) term) ? Typically γ=0.9 if constraint is violated and γ=1 otherwise ? Doesn’t affect convergence, since specific form of constraint is only crucial when x k is close to x* T 1 () ()0 1,, kk k jjj ggjmγ?+≤=xS x … 30 ? Massachusetts Institute of Technology - Prof. de Weck and Prof. Willcox Engineering Systems Division and Dept. of Aeronautics and Astronautics Sequential Quadratic Programming ? Widely used in engineering applications e.g. NLOpt ? Considered to be the best gradient-based algorithm ? Strong theoretical basis 31 ? Massachusetts Institute of Technology - Prof. de Weck and Prof. Willcox Engineering Systems Division and Dept. of Aeronautics and Astronautics Mixed Integer Programming ? Rounding – unlikely to achieve optimum discrete solution – rounded solution may be infeasible ? Branch and Bound – create tree with several optimization branches – continuous solution is a lower bound on the mixed solution – usually solve a large number of optimization problems – expensive 32 ? Massachusetts Institute of Technology - Prof. de Weck and Prof. Willcox Engineering Systems Division and Dept. of Aeronautics and Astronautics Lecture Summary We have seen the following nonlinear techniques: ? Penalty and Barrier Methods ? Sequential Quadratic Programming It is important to understand ? when it is appropriate to use these methods ? the basics of how the method works ? why the method might fail ? what your results mean (numerically vs. physically) 33 ? Massachusetts Institute of Technology - Prof. de Weck and Prof. Willcox Engineering Systems Division and Dept. of Aeronautics and Astronautics References Practical Optimization, P.E. Gill, W. Murray and M.H. Wright, Academic Press, 1986. Numerical Optimization Techniques for Engineering Design, G.N. Vanderplaats, Vanderplaats R&D, 1999. Optimal Design in Multidisciplinary Systems, AIAA Professional Development Short Course Notes, September 2002.