Exercise session 3

1. Exercise session 3 Assoc. Prof. Dr.Pelin Gundes

2. Numerical techniques for constrained optimization There are two outcomes that the algorithms seek to accomplish: • The first is to ensure that the design is feasible (satisfies all constraints) • The second is that it is optimal (it satisfies the Kuhn-Tucker conditions) • While the focus is on determining the solution, in times of difficulty, it is essential to remember that feasibility is more important than optimality. • Two distinct approaches will be used to handle the constrained optimization problem. • The first approach is termed the indirect approach and solves the problem by transforming it into an unconstrained problem. • The second approach is to handle the constraints without transformation- the direct approach. • Two indirect methods are presented: The Exterior Penalty Function Methods and the Augmented Lagrange Multiplier Method. • The direct approach handles the constraints and the objective together without any transformation. Four methods are presented: The Sequential Linear Programming, Sequential Quadratic Programming, Generalized Reduced Gradient Method and Sequential Gradient Restoration Algorithm.

3. Direct methods for constrained optimization • The direct methods include both the objective and the constraints to search for the optimal solution. • Most of them are based on linearization of the functions about the current design point. Linearization is based on expansion of the function about the current variable values using the Taylor series. • Linearization:The Taylor series for a two variable expanded function f(x,y) quadratically about the current point (xp ,yp) is expressed as:

4. Direct methods for constrained optimization • If the displacements are organized as a column vector, the expansion can be expressed in a condensed manner as: • For n variables,with Xp the current point and ΔX the displacement vector, • The above equation can be written in terms of the difference in function values as:

5. Direct methods for constrained optimization • Where f=fTΔX is termed the first variation, 2f is the second variation and is given by the second term in the above equation. • In linearization of the function f about the current value of design Xp, only the first variation is used. The neighboring value of the function can be expressed as: • All of the functions in the problem can be linearized similarly. It is essential to understand the difference between f (X) and

6. Direct methods for constrained optimization • This is illustrated in the figure using the objective function expanded about the current design x1=3 and x2=2. Thecurved linesf (X) are thecontours of the objective function below. The straight lines are the contours of the linearized function (lines of constant value of the function) obtained through the following:

7. Direct methods for constrained optimization • The quadratic expansion of the function can be obtained by using • The expanded curves will be nonlinear (quadratic). They appear as ellipses in the figure where several contours of the objective function are expanded quadratically about x1=3 and x2=2.

8. Direct methods for constrained optimization Four direct methods that will be discussed are: • Sequential Linear Programming where the solution is obtained by successively solving the corresponding linearized optimization problem. • Sequential Quadratic Programming uses the quadratic expansion for the objective function. Like the SUMT, the current solution provides the starting values for the next iteration. • Generalized Reduced Gradient Method develops a sophisticated search direction and follows it nwith an elaborate one-dimensional process. • Sequential Gradient Restoration Algorithm uses a two cycle approach working on feasibility and optimality alternately to find the optimum.

9. Direct methods for constrained optimization

10. Direct methods for constrained optimization

11. Direct methods for constrained optimization

12. Direct methods for constrained optimization

13. Direct methods for constrained optimization % % Optimization with MATAB; Dr P.Venkataramana= % Chapter 7 Section 7.3.1 % Sequential Linear Programming – % Generating plots for a chosen design vector % Example 7.1 % Symbolic Calculations and plotting % Two variable problem only % format compact syms x1 x2 f g h syms gradf1 gradf2 gradh1 gradh2 gradg1 gradg2 f = x1^4 - 2*x1*x1*x2 + x1*x1 + x1*x2*x2 - 2*x1 + 4; h = x1*x1 + x2*x2 - 2; g = 0.25*x1*x1 +0.75*x2*x2 -1; gradf1 = diff(f,x1); gradf2 = diff(f,x2); gradh1 = diff(h,x1); gradh2 = diff(h,x2); gradg1 = diff(g,x1); gradg2 = diff(g,x2);

14. Direct methods for constrained optimization string1 = ['Input the design vector: ']; xb = input(string1); fprintf('Initial design vector\n'),disp(xb) f1 = subs(f,{x1,x2},{xb(1),xb(2)}) g1 = subs(g,{x1,x2},{xb(1),xb(2)}) h1 = subs(h,{x1,x2},{xb(1),xb(2)}) gf1 = double(subs(gradf1,{x1,x2},{xb(1),xb(2)})); gf2 = double(subs(gradf2,{x1,x2},{xb(1),xb(2)})); gh1 = double(subs(gradh1,{x1,x2},{xb(1),xb(2)})); gh2 = double(subs(gradh2,{x1,x2},{xb(1),xb(2)})); gg1 = double(subs(gradg1,{x1,x2},{xb(1),xb(2)})); gg2 = double(subs(gradg2,{x1,x2},{xb(1),xb(2)}));

15. Direct methods for constrained optimization % choose values for rh and rg x11 = -5:.2:5; x22 = -5:.2:5; x1len = length(x11); x2len = length(x22); for i = 1:x1len; for j = 1:x2len; xbv = [x11(i) x22(j)]; fnew(j,i) = f1 + [gf1 gf2]*xbv'; hnew(j,i) = h1 + [gh1 gh2]*xbv'; gnew(j,i) = g1 + [gg1 gg2]*xbv'; end end

16. Direct methods for constrained optimization minf = min(min(fnew)); maxf = max(max(fnew)); mm = (maxf - minf)/5.0; mvect=[(minf+mm) (minf + 1.5*mm) (minf + 2*mm) (minf+ 2.5*mm) (minf + 3.0*mm) (minf +4*mm) (minf + 4.5*mm)]; [c1,fc] = contour(x11,x22,fnew,mvect,'b'); clabel(c1); set(fc,'LineWidth',2) Grid xlabel('\delta x_1') ylabel('\delta x_2') title('Example 7.1: Sequential Linear Programming') hold on

17. Direct methods for constrained optimization [c2,hc]=contour(x11,x22,hnew,[0,0],'g'); set(hc,'LineWidth',2,'LineStyle','--') grid [c3,gc]=contour(x11,x22,gnew,[0,0],'r'); set(gc,'LineWidth',2,'LineStyle',':') contour(x11,x22,gnew,[1,1],'k') grid hold off

18. Indirect methods for constrained optimization • These methods were developed to take advantage of codes that solve unconstrained optimization problems. They are also referred to as Sequential Unconstrained Minimization Techniques (SUMT). • The idea behind the approach is to repeatedly call the unconstrained optimization algorithm using the solution of the previous iteration. The unconstrained algorithm itself executes many iterations. This would require a robust unconstrained minimizer to handle a large class of problems. The BFGS method is robust and impressive over a large class of problems.

19. Indirect methods for constrained optimization • A preprocessing task involves transforming the constrained problem into an unconstrained problem. This is largely accomplished by augmenting the objective function with additional functions that reflect the violation of the constraints very significantly. These functions are referred to as the penalty functions. • There was significant activity in this area at one time which led to entire families of different types of penalty function methods.The first of these was the Exterior Penalty Function Method. The first version of ANSYS that incorporated optimization in its finite element program relied on EFP. The EFP had several short comings. • To address those, the interior penalty function methods were developed leading to the Variable Penalty Function methods. • In this chapter, only the EPF is addressed largely due to academic interest. In view of the excellent performance of the direct methods, these methods will probably not be used today for continuous problems. They are once again important in global optimization techniques for constrained problems.

20. Indirect methods for constrained optimization • The second method presented in this section, the AugmentedLagrange Method (ALM) is the best of the Sequential Unconstrained Minimization Techniques. Its exceedingly simple implementaion, its quality of solution, and its ability to generate information on the Lagrange multipliers allow it to seriously challenge the direct techniques.

21. Exterior Penalty Function Method • The transformation of the optimization problem to an unconstrained problem is made possible through a penalty function formulation. The transformed unconstrained problem is: where P(X, rh, rg) is the penalty function, rh and rg are penalty constants (also called multipliers).

22. Exterior Penalty Function Method • The penalty function is expressed as: In the above equation, if the equality constraints are not zero, their value gets squared and multiplied by the penalty parameter and then gets added to the objective function. If the inequality constraint is in violation, it too gets squared and added to the objective function after being amplified by the penalty multipliers. In a sense if the constraints are not satisfied, then they are penalized, hence the function’s name. • It can be shown that the transformed unconstrained problem solves the original constrained problem as the multipliers rh and rg approach . • In computer implementation,this limit is replaced by a large value instead of . Another facet of the computer implementation of this method is that a large value of the multipliers at the first iteration is bound to create numerical difficulties.

23. Exterior Penalty Function Method • These multipliers are started with small values and updated geometrically with each iteration. The unconstrained technique, for example DFP, will solve the equation below for a known value of the multipliers. • The solution returned from DFP, will solve the above equation for a known value of the multipliers.The solution returned from the DFP can be considered as a function of the multiplier and can be thought of as: X*=X*(rh ,rg ) • The Sequential Unconstrained Minimization Techniques (SUMT) iteration involves updating the multipliers and the initial design vector and calling the unconstrained minimizer again.

24. Exterior Penalty Function Method

25. Exterior Penalty Function Method

26. Exterior Penalty Function Method • The EPF is very sensitive to the starting value of the multipliers and to the scaling factors as well. Different problems respond favorably to different values of the multipliers. • It is recommended that the initial values of the multipliers be chosen as the ratio of the objective function to the corresponding term in the penalty function at the initial design. This ensures that both the objective function and the constraints are equally important in determining the changes in the design for the succeeding iteration. • One reason for the term Exterior Penalty is that at the end of each SUMT iteration, the design will be infeasible (until the solution is obtained). This implies that the method determines design values that are approaching the feasible region from the outside. This is a serious drawback if the method fails prematurely, as it will often do. The information generated so far is valueless as the designs were never feasible.

27. Exterior Penalty Function Method • As seen in the example below, the EPF severely increases the nonlinearity of the problem creating conditions for the method to fail. • It is expected that the increase in nonlinearity is balanced by a closer starting value for the design, as each SUMT iteration starts closer to the solution than the previous one. • In the following slides, the EFP method is applied to an example through a series of calculations rather than through the translation of the algorithm into MATLAB code. • There are a couple of changes with respect to algorithm A7.2. To resolve the penalty function with respect to the inequality constraint, the constraint is assumed to always be in violation so that the return from the max function is the constraint function itself. This will drive the inequality constraint to be active which we know to be true for this example. • Numerical implementation as outlined in the algorithm should allow the determination of inactive constraints. • Instead of the numerical implementation of the unconstrained problem, an analytical solution is determined using MATLAB symbolic computation.

28. Exterior Penalty Function Method

29. Exterior Penalty Function Method • The figure below is the contour plot of the original graphical solution for the example.

30. Exterior Penalty Function Method • The figure below is the contour plot of the transformed unconstrained function for values of rh=1 and rg=1 . The increase in nonlinearity is readily apparent.

31. Exterior Penalty Function Method • The figure below is the plot of rh=5 and rg=5 . This and the previous figure suggest several points that satisfy first order conditions. Their closeness makes it difficult for any numerical technique to find the optimum. It is clear that the EPF severely increases the nonlinearity of the problem.

32. Exterior Penalty Function Method • The code below uses symbolic manipulation for the exterior penalty function method. Since the code uses symbolic manipulation, actually drawing the plot takes some time. Evaluating the data numerically will make a big difference. • Before applying the code, the corresponding unconstrained function for this problem is constructed: • The code is an m file that will calculate the solution and the values of the function for a predetermined set of values of the penalty multipliers.

33. Exterior Penalty Function Method • The code requires two inputs from the user at different stages. The first input is the values for the multipliers for which the solution should be obtained. • The list of solution (there are nine for this problem) is displayed in the Command Window. • The user finds the solution that satisfies the side constraints which must be entered at the prompt. (note: it must be entered as a vector). The program then prints out the values of the various functions involved in the example).

34. Exterior Penalty Function Method % % Optimization with MATAB; Dr P.Venkataramana= % Chapter 7 Section 7.2.1 % External Penalty Function Method % Example 7.1 % Symbolic Calculations and plotting format compact syms x1 x2 rh rg f g h F grad1 grad2 f = x1^4 - 2*x1*x1*x2 + x1*x1 + x1*x2*x2 - 2*x1 + 4; h = x1*x1 + x2*x2 - 2; g = 0.25*x1*x1 +0.75*x2*x2 -1; %F = f + rh*h*h + rg*g*g; %grad1 = diff(F,x1); %grad2 = diff(F,x2); % choose values for rh and rg rh = 5; rg = 5;

35. Exterior Penalty Function Method x11 = -2:.2:5; x22 = -2:.2:5; x1len = length(x11); x2len = length(x22); for i = 1:x1len; for j = 1:x2len; gval = subs(g,{x1 x2},{x11(i) x22(j)}); if gval < 0 gval = 0; end hval = subs(h,{x1 x2},{x11(i) x22(j)}); Fval(j,i) = subs(f,{x1 x2},{x11(i) x22(j)}) ... + rg*gval*gval + rh*hval*hval; end end c1 = contour(x11,x22,Fval,[3.1 4 5 6 10 20 50 100 200 500]); clabel(c1); grid xlabel('x_1') ylabel('x_2') strng = strcat('Example 7.1: ','r_h = ', num2str(rh),'r_g = ',num2str(rg)); title(strng)

36. Exterior Penalty Function Method • The following is posted from the command window for both of the penalty multipliers set to 25.

37. Exterior Penalty Function Method • In the above run, the equality and the inequality constraints are not satisfied. Another m file for implementing the Exterior Penalty Function method is pasted below: % % Optimization with MATLAB; Dr P.Venkataramana= % Chapter 7 Section 7.2.1 % External Penalty Function Method % Symbolic Calculations (partial) % Example 7.1 format compact % define the functions syms x1 x2 rh rg f g h F grad1 grad2 f = x1^4 - 2*x1*x1*x2 + x1*x1 + x1*x2*x2 - 2*x1 + 4; h = x1*x1 + x2*x2 - 2; g = 0.25*x1*x1 +0.75*x2*x2 -1;

38. Exterior Penalty Function Method % the unconstrained function F = f + rh*h*h + rg*g*g; % gradients of the unconstrained function grad1 = diff(F,x1); grad2 = diff(F,x2); % choose values for rh and rg fprintf('\n') rh = input('enter value for rh [default = 1] : '); if isempty(rh) rh = 1; end fprintf('\n') rg = input('enter value for rg [default = 1] : '); if isempty(rg) rg = 1; end

39. Exterior Penalty Function Method % solve for x1 and x2 sol = solve(subs(grad1), subs(grad2)); %display all the solutions to choose 1 [double(sol.x1) double(sol.x2)] % enter the value for design vector scanned from command window string1 = ['Input the design vector chosen for evaluation.\n'] ; xs = input(string1); fprintf('\nThe design vector [ %10.4f %10.4f ]',xs); fv = subs(f,{x1,x2},{xs(1),xs(2)}); hv = subs(h,{x1,x2},{xs(1),xs(2)}); gv = subs(g,{x1,x2},{xs(1),xs(2)}); fprintf('\nobjective function: '), disp(fv) fprintf('\nequality constraint: '), disp(hv) fprintf('\ninequality constraint: '), disp(gv)

40. Exterior Penalty Function Method • In the above, the scaling factor is 5 for both multipliers. A glance at the Table below clearly illustrates the characteristics of the EPF method.

41. Exterior Penalty Function Method As the values of the multipliers increase: • The design approaches the optimal value. • The constraint violations decrease. • The solution is being approached from outside the feasible region.

42. Augmented Lagrangian Method • This is the most robust of the penalty function methods. • More importantly, it also provides information on the Lagrange multipliers at the solution. This is achieved by not solving for the multipliers but merely updating them during successive SUMT iterations. • It overcomes many of the difficulties associated with the penalty function formulation without any significant overhead.

43. Augmented Lagrangian Method Transformation of the Unconstrained Problem: The general optimization problem restated below EQUATION 7.1 TO 7.4 is transformed as in the method of Lagrange multipliers EQUATION 7.23 TO 7.24

44. Augmented Lagrangian Method • Here  is the multiplier vector tied to the equality constraints, is the multiplier vector associated with the inequality constraints, and rhand rg are the penalty multipliers used similar to the EPF method. • F is solved as an unconstrained function for predetermined values of , , rhand rg . Therefore, the solution for each SUMT iteration is X*=X*(, , rh, rg) • At the end of the SUMT iteration, the values of the multipliers and penalty constants are updated. The latter are usually geometrically scaled but unlike EPF do not have to be driven to  for convergence.

45. Augmented Lagrangian Method ALM.m code