Lus F. D. Brs Alvaro F. M. Azevedo Faculty of Engineering - PowerPoint PPT Presentation

OBJECT ORIENTED IMPLEMENTATION OF A SECOND-ORDER OPTIMIZATION METHOD Luís F. D. Brás Alvaro F. M. Azevedo Faculty of Engineering University of Porto PORTUGAL 28-30 May 2001 OPTI 2001 Bologna - Italy 1

OPTIMIZATION APPROACH • Nonlinear program • Second-order approximation • Integrated formulation • All the problem variables are present in the nonlinear program • No sensitivity analysis 28-30 May 2001 OPTI 2001 Bologna - Italy 2

OPTIMIZATION SOFTWARE • General purpose code • Lagrange-Newton method • Symbolic manipulation of all the functions • Exact 1 st and 2 nd derivatives • Object oriented approach • Language: C++ 28-30 May 2001 OPTI 2001 Bologna - Italy 3

OBJECT ORIENTED PROGRAMMING • What we gain ♦ Higher abstraction level ♦ Encapsulation of lower level complexities ♦ Code maintenance and reuse is facilitated • What we lose ♦ Performance ♦ Straightforward coding 28-30 May 2001 OPTI 2001 Bologna - Italy 4

OBJECT ORIENTED FEATURES • Classes • Function and operator overloading • Inheritance • Polymorphism • Templates • Exception handling 28-30 May 2001 OPTI 2001 Bologna - Italy 5

NONLINEAR PROGRAMMING M inim ize f ( ) x ~ subject to → = ≤ 0 2 g ( ) + x s 0 g x ( ) i i ~ ~ ~ ~ = 0 h x ( ) ~ ~ ~ • Variables / functions real and continuous • Generic functions can be treated 28-30 May 2001 OPTI 2001 Bologna - Italy 6

LAGRANGIAN ( ) ( ) ( ) ( ) p m   ∑ ∑ = + λ +  + λ g 2 h L X f x g x s h x    k k k k k ~ ~ ~ ~ = = k 1 k 1 VARIABLES ( ) = λ λ g h X x , s , , ~ ~ ~ ~ ~ SOLUTION • Stationary point of the Lagrangian 28-30 May 2001 OPTI 2001 Bologna - Italy 7

SYSTEM OF NONLINEAR EQUATIONS ( ) = 1,..., λ = g i m 2 0 s i i ( ) = 1,..., + = 2 i m 0 g s ( ) i i ∇ = ⇒ L X 0 ∂ ∂ ∂ p m f g h ( ) ∑ ∑ = 1,..., + λ + λ = ~ ~ g h k k i n 0 ∂ ∂ ∂ k k x x x = = k 1 k 1 i i i ( ) = 1,..., h i = 0 i p • The solution of the system is a KKT solution when λ g ≥ 0 ~ 28-30 May 2001 OPTI 2001 Bologna - Italy 8

LAGRANGE-NEWTON METHOD • The system of nonlinear equations ∇ = L X ( ) 0 ~ ~ is solved by the Newton method • In each iteration the following system of linear equations has to be solved ( ) ( ) − − ∆ + ∇ = q 1 q q 1 H X X L X 0 ~ ~ ~ ~ ~ • H is the Hessian of the Lagrangian • Second derivatives of all the functions are required 28-30 May 2001 OPTI 2001 Bologna - Italy 9

SYSTEM OF LINEAR EQUATIONS • Gaussian elimination ♦ adapted to the sparsity pattern of the Hessian matrix • Conjugate gradients ♦ diagonal preconditioning ♦ adapted to an indefinite Hessian matrix 28-30 May 2001 OPTI 2001 Bologna - Italy 10

LINE SEARCH − 1 = + α ∆ q q q X X X ~ ~ ~ • The value of α minimizes the error in ∆ X q direction ~ ♦ the value of α is often close to one ♦ faster convergence ♦ process may fail • The value of α is made considerably smaller (e.g. α = 0.1) ♦ stable convergence ♦ more iterations - slower 28-30 May 2001 OPTI 2001 Bologna - Italy 11

AUTOMATIC DIFFERENTIATION • Expression evaluation • Partial derivative calculation (first, second, ...) • Each function is parsed and stored as a tree of tokens (constants, variables and operators) • Automatic differentiation is based on Rall numbers 28-30 May 2001 OPTI 2001 Bologna - Italy 12

RALL NUMBERS • A Rall number is a class that encapsulates the numerical value of the function, its gradient vector and its Hessian matrix • All the operators are overloaded in order to apply the differentiation rules • With Rall numbers automatic differentiation can be efficiently performed 28-30 May 2001 OPTI 2001 Bologna - Italy 13

RALL NUMBERS Example: ( ) ( ) f x 1 , x g x 1 , x Functions and 2 2 Derivatives of the product: ∂ ∂ ∂ ( ) f g = + f g g f ∂ ∂ ∂ x x x 1 1 1 ∂ ∂ ∂ ∂ ∂ ∂ ∂ 2 2 2 ( ) f f g f g g = + + + f g g f ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ x x x x x x x x x x 1 2 1 2 2 1 1 2 1 2 28-30 May 2001 OPTI 2001 Bologna - Italy 14

RALL NUMBERS class CRall { double x; // Operand value double v[2]; // df/dx1, df/dx2 double m[2][2]; // d2f/dxi dxj public: CRall CRall::operator* (const CRall & g) const { CRall t; t.x = x * g.x; t.v[0] = v[0]*g.x + x*g.v[0]; t.v[1] = v[1]*g.x + x*g.v[1]; t.m[0][0] = m[0][0]*g.x+v[0]*g.v[0]+v[0]*g.v[0]+x*g.m[0][0]; t.m[0][1] = m[0][1]*g.x+v[0]*g.v[1]+v[1]*g.v[0]+x*g.m[0][1]; t.m[1][0] = m[1][0]*g.x+v[1]*g.v[0]+v[0]*g.v[1]+x*g.m[1][0]; t.m[1][1] = m[1][1]*g.x+v[1]*g.v[1]+v[1]*g.v[1]+x*g.m[1][1]; return t; } }; 28-30 May 2001 OPTI 2001 Bologna - Italy 15

RALL NUMBERS x = constant value; v = [0,0]; m = [[0,0],[0,0]] x = value of x 1 ; v = [1,0]; m = [[0,0],[0,0]] x = value of x 2 ; v = [0,1]; m = [[0,0],[0,0]] 28-30 May 2001 OPTI 2001 Bologna - Italy 16

EXPRESSION PARSER • A binary tree is constructed according to the operator precedence • Each tree node is a Rall number • A symbol table is initialized with the values of the variables and constants • The tree traversal causes an evaluation of the function, gradient and Hessian 28-30 May 2001 OPTI 2001 Bologna - Italy 17

EXPRESSION PARSER ( ) = + − 2 • Example: f x , x , x ( x 8 x ) / ( 6 x ) 1 2 3 1 2 3 x 2 x 3 8 2 6 x 1 * ^ − + / 28-30 May 2001 OPTI 2001 Bologna - Italy 18

EXPRESSION PARSER Exp_Obj Binary_Obj Base_Expression Mul_Obj Pow_Obj CExpression Binary_Add_Obj Close_Paren_Obj Binary_Sub_Obj Literal_Obj Div_Obj Open_Paren_Obj Start_Obj Unary_Obj Stop_Obj Unary_Log10 Unary_Minus Variable_Obj Unary_Loge Unary_Sin Unary_Cos Unary_Sqr ExpStack Unary_Exp Unary_Sqrt Symbol_Table 28-30 May 2001 OPTI 2001 Bologna - Italy 19

SCALING x = c x • Variable substitution: i i g = • Constraint normalization: c g i i Min . 2000 x Min . y 1 1 subject to subject to − + + = − + + = 2 2 0 . 640 0 . 256 0 . 384 0 x 200 x 0 y y 1 3 1 3 − + = − + = 2 2 x 0 . 2 x 0 0 . 447 y 0 . 894 0 . 447 y 0 2 4 2 4 − + = − + = 0 . 707 y y 0 . 707 0 10 x x 500 0 1 2 1 2 28-30 May 2001 OPTI 2001 Bologna - Italy 20

NUMERICAL EXAMPLE x 2 F = (f 1 ,f 2 ) 3 A A (cross-sectional area) 100 cm r = 2 0 0 F kN = 1 2 f 8 f x 1 2 1 σ m ax = 10 0 2 kN cm x x • Variables: A , x • Svanberg’s solution confirmed 28-30 May 2001 OPTI 2001 Bologna - Italy 21

NONLINEAR PROGRAM ( ) = + 2 Min . w x , x C x 1 x 1 2 1 1 2 subject to   ( ) 8 1   σ = + + ≤ 2 , 1 1 x x C x   1 1 2 2 2   x x x 1 1 2   ( ) 8 1   σ = + − ≤ 2 x , x C 1 x 1   2 1 2 2 2   x x x 1 1 2 ≤ ≤ ≤ ≤ 0 . 2 x 4 . 0 ; 0 . 1 x 1 . 6 1 2 28-30 May 2001 OPTI 2001 Bologna - Italy 22

DATA FILE # Main title # Minimum x1 Shape optimization of a two bar truss -x1 + 0.2; # N. of eq. constr.; N. of ineq. constr. # Maximum x1 0 6 x1 - 4.0; # Objective Function # Minimum x2 C1 * x1 * sqrt(1+x2^2); -x2 + 0.1; # Allowable stress - bar 1 # Maximum x2 C2 * sqrt(1+x2^2) * (8/x1+1/x1/x2) - 1; x2 - 1.6; # Allowable stress - bar 2 # N. of variables C2 * sqrt(1+x2^2) * (8/x1-1/x1/x2) - 1; 4 SUBSTITUTED, 1.000, C1; SUBSTITUTED, 0.124, C2; INDEPENDENT, 1.5, x1; INDEPENDENT, 0.5, x2; 28-30 May 2001 OPTI 2001 Bologna - Italy 23

38 bar truss - member sizing Stress and displacement constraints 28-30 May 2001 OPTI 2001 Bologna - Italy 24

CONCLUSIONS • Code maintenance • Efficiency and accuracy in the evaluation of derivatives • Easy inclusion of alternative numerical techniques • Not efficient in the OO manipulation of the Hessian matrix • Friendly user interface is still required 28-30 May 2001 OPTI 2001 Bologna - Italy 25