Introduction to Optimizatio A Brief Review of Linear Algebra

z A brief review of Linear Algebra z Another LP formulation, and z Math Programming and Radiation Therapy 1 Review of Linear Algebra z Some elementary facts about vectors and matrices. z The Gauss-Jordan method for solving systems of equations. z Bases and basic solutions and pivoting. Note: you are not expected to know 18.06, but will be expected to know about solving systems of equations. 2 Elementary Facts about Vectors v = ⎡⎣v1 v2 v3 v4 ⎤⎦ is called a row vector. ⎡ v1 ⎤ ⎢v ⎥ t ⎢ 2⎥ v = The transpose of v is a column vector. ⎢ v3 ⎥ ⎢ ⎥ ⎣ v4 ⎦ w = [ w1 w2 w3 w4 ] is another row vector. The inner product of vectors v and w is given by: v D w = vw = v1w1 + v 2 w2 + v3 w3 + v4 w4 t 3 Matrix Multiplication A=(a ) ij B =(b ) ij C =(c ) = A× B ij Suppose that A has n columns and B has n rows. cij = ∑k=1aikbkj n 4 Multiplying Matrices Let C = (cij) = A × B. Then cij is the inner product of row i of A and column j of B. ⎡ a11 A = ⎢⎢ a21 ⎢⎣ a31 a12 a22 a32 a13 ⎤ ⎡ b11 a23 ⎥⎥ B = ⎢⎢ b21 ⎢⎣ b31 a33 ⎥⎦ b12 b22 b32 b13 ⎤ b23 ⎥⎥ b33 ⎥⎦ For example, what is c23? c23 = a21b13 + a22 b23 + a23 b33 5 Multiplying Matrices 2 Let C = (cij) = A × B. Then each column of C is obtained by adding multiples of columns of A. ⎡ 1 2 3 ⎤ ⎡100 ⎤ ⎡ ⎢ 4 5 6 ⎥ × ⎢ 10 ⎥ = ⎢ ⎢ ⎥ ⎢ ⎥ ⎢ ⎢⎣ 7 8 9 ⎥⎦ ⎢⎣ 1 ⎥⎦ ⎢⎣ ⎡1⎤ ⎡ 2⎤ ⎡ 3⎤ = 100 ⎢ 4 ⎥ + 10 ⎢ 5 ⎥ + ⎢ 6 ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢⎣ 7 ⎥⎦ ⎢⎣ 8 ⎥⎦ ⎢⎣ 9 ⎥⎦ ⎤ ⎥ ⎥ ⎥⎦ Take 30 seconds and obtain the answer. This is not a coincidence. Matrix multiplication was defined so that this would be true. 6 Multiplying Matrices 3 Let C = (cij) = A × B. Then each row of C is obtained by adding multiples of rows of A. ⎡1 2 3⎤ ⎢4 5 6⎥ = × 100 10 1 [ ] ⎢ ] ⎥ [ ⎢⎣ 7 8 9 ⎥⎦ Take 30 seconds and obtain the answer. = 100 [1 2 3] + 10 [ 4 5 6] + [ 7 8 9] Matrix Premultiplication ⇔ Adding multiples of rows. 7 Elementary Facts about Solving Equations Solve for Ax = b, where ⎡ x1 ⎤ ⎡1 2 4 ⎤ x = ⎢ x ⎥ A=⎢ 2⎥ ⎢ ⎥ ⎣ 2 1 −1⎦ ⎢⎣ x3 ⎥⎦ 2×3 3×1 x1 + 2 x2 + 4 x3 = 0 2 x1 + x2 − x3 = 6 Find a linear combination of the columns of A that equals b. ⎡1⎤ ⎡ 2⎤ ⎡4⎤ ⎡ 0⎤ ⎢ 2 ⎥ x1 + ⎢ 1 ⎥ x2 + ⎢ −1⎥ x3 = ⎢ 6 ⎥ ⎣ ⎦ ⎣ ⎦ ⎣ ⎦ ⎣ ⎦ ⎡0⎤ b=⎢ ⎥ ⎣6⎦ 2×1 Solving a System of Equations x1 x2 x3 x4 1 2 -1 2 1 1 4 -1 2 1 -1 2 = = = 0 6 -3 To solve a system of equations, use Gauss-Jordan elimination. 9 To solve the system of equations: x1 x2 x3 x4 1 2 4 1 = 0 2 1 -1 -1 = 6 -1 1 2 2 = -3 G-J Animation 10 The fundamental operation: pivoting a special case of an elementary row operation (ero) x1 x2 x3 x4 a11 a12 a13 a14 = b1 a21 a22 a23 a24 = b2 a31 a32 a33 a34 = b3 Pivot on a23 11 Pivot on a23 x1 x2 x3 x4 a11 a12 a13 a14 = b1 a21a/a 21 23 a22 a123 a24 = b2 a31 a32 a33 a34 = b3 Take two minutes to fill in the other coefficients for constraint 2 and constraint 1 after the pivot. 12 Jordan Canonical Form for an m x n matrix x1 x2 x3 x4 1 0 0 -1 = 2 0 1 0 1 = 1 0 0 1 0 = -1 There are m columns that have been transformed into unit vectors, one for each row. The variables in these columns are called “basic.” The remaining variable(s) are called non-basic. 13 Jordan Canonical Form for an m x n matrix x1 x2 x3 x4 1 0 0 -1 = 2 0 1 0 1 = 1 0 0 1 0 = -1 In the “basic solution”, the non-basic variables are set to 0. The basic solution is: x1 = 2, x2 = 1, x3 = -1 x4 = 0 14 There is an easily determined solution for every choice of non-basic variables. Suppose we set x4 = 2. x1 x2 x3 x24 1 0 0 -1 = 2 0 1 0 1 = 1 0 0 1 0 = -1 Equation 1 becomes: x1 – 2 = 2. So x1 = 4. If we set x4 = 2, we get x1 = 4, x2 = -1, x3 = -1 x4 = 2. 15 Another Jordan Canonical Form x1 x2 x3 x4 0 1 -2 0 -1 1 = 4 1 1 1 3 0 1 0 = -2 3 0 0 -3 1 0 = -1 3 Take 60 seconds and answer the following questions. What are the “basic” variables? What is the basic solution? What is the solution if the non-basic variable is set to 2? 16 The pigskin problem (from Practical Management Science) z Pigskin company makes footballs z All data below is for 1000s of footballs z Forecast demand for next 6 months – 10, 15, 30, 35, 25 and 10 z Current inventory of footballs: 5 z Determine the production levels and inventory levels over the next six months – meet demand at minimum cost – satisfy constraints 17 The pigskin problem (continued) z Max Production capacity: 30 per month z Max Storage capacity: 10 per month z Production Cost per football for next 6 months: – $12.50, $12.55, $12.70, $12.80, $12.85, $12.95 z Holding cost: $.60 per football per month z With your partner: write an LP to describe the problem, but only write variables and constraints for the next 3 months. 18 On the formulation z Choose decision variables. – Let xj = the number of footballs produced in month j (in 1000s) – Let yj = the number of footballs held in inventory from month j to month j + 1. (in 1000s) – y0 = 5 z Then write the constraints and the objective. 19 The formulation Minimize subject to y0 = 5 0 <= yj <= 10 for j = 1, 2, 3 0 <= xj <= 30 for j = 1, 2, 3 Pigskin Spreadsheet 20 Math Programming and Radiation Therapy z An important application area for optimization z 21 Math Programming and Radiation Therapy z High doses of radiation (energy/unit mass) can kill cells and/or prevent them from growing and dividing – True for cancer cells and normal cells z Radiation is attractive because the repair mechanisms for cancer cells is less efficient than for normal cells 22 Radiation Therapy Overview z Recent advances in radiation therapy now make it possible to – map the cancerous region in greater detail – aim a larger number of different beamlets with greater specificity z This has spawned the new field of tomotherapy z “Optimizing the Delivery of Radiation Therapy to Cancer patients,” 23 Conventional Radiotherapy Relative Intensity of Dose Delivered 24 Conventional Radiotherapy Relative Intensity of Dose Delivered 25 Conventional Radiotherapy z In conventional radiotherapy – 3 to 7 beams of radiation – radiation oncologist and physicist work together to determine a set of beam angles and beam intensities – determined by manual “trial-and-error” process 26 Goal: maximize the dose to the tumor while minimizing dose to the critical area Critical Area Tumor area With a small number of beams, it is difficult to achieve these goals. 27 Recent Advances z More accurate map of tumor area – CT -- Computed Tomography – MRI -- Magnetic Resonance Imaging z More accurate delivery of radiation – IMRT: Intensity Modulated Radiation Therapy – Tomotherapy 28 Tomotherapy: a diagram 29 Radiation Therapy: Problem Statement z For a given tumor and given critical areas z For a given set of possible beamlet origins and angles z Determine the weight of each beamlet such that: – dosage over the tumor area will be at least a target level γL . – dosage over the critical area will be at most a target level γU. 30 Display of radiation levels 31 Linear Programming Model z First, discretize the space – Divide up region into a 2D (or 3D) grid of pixels 32 More on the LP z Create the beamlet data for each of p = 1, ..., n possible beamlets. z Dp is the matrix of unit doses delivered by beam p. p ij D = unit dose delivered to pixel (i, j) by beamlet p 33 Linear Program z Decision variables w = (w1, ..., wp) z wp = intensity weight assigned to beamlet p for p = 1 to n; z Dij = dosage delivered to pixel (i, j) Dij = ∑ p=1 D w p n p ij 34 An LP model minimize ∑ D ij (i, j) took 4 minutes to solve. Dij = ∑ p=1 D w p n p ij Dij ≥ γ L for ( i , j ) ∈ T Dij ≤ γ U for ( i , j ) ∈ C wp ≥ 0 for all p In an example reported in the paper, there were more than 63,000 variables, and more than 94,000 constraints (excluding upper/lower bounds) 35 What happens if the model is infeasible? minimize ∑ D y ((ii,, jj)) ijij Dij = ∑ p=1 D w p n p ij Allow the constraint for pixil (i,j) to be violated by an amount yij, and then minimize the violations. ( i , j()i ∈ DD + ≥yijγ L≥ γ L for for , jT ) ∈T ij ij DD ( i , j()i∈ , jC )∈C ij ij− ≤yγ ij U≤ γ U for for wp ≥ 0 yij ≥ 0 for all p for all ( i , j ) 36 An even better model minimize D y y ∑ ( ) ∑ 2 ijij ( i ,((jii),, jj)) ij Dij = ∑ p=1 D w p n minimize the sum of squared violations. p ij ( i , j()i ∈ DD + ≥yijγ L≥ γ L for for , jT ) ∈T ij ij DD ( i , j()i∈ , jC )∈C ij ij− ≤yγ ij U≤ γ U for for wp ≥ 0 yij ≥ 0 Least squares for all p for all ( i , j ) This is a nonlinear program (NLP). This one can be solved efficiently. 37 Optimal Solution for the LP 38 An Optimal Solution to an NLP 39 Is that the end of the story on modeling? z Other issues: – Delivering the radiation doses quicker by setting the multi-leaf collimators optimally. – trying to keep average radiation levels low over part of the critical region. – how do we trade off the doses to the critical area with doses to the tumor? 40 The 2003 Edelman Prize for OR/MS Practice Winner: Canadian Pacific Railway "Perfecting the Scheduled Railroad: Model-Driven Operating Plan Development." – Previous operation: Waits for enough cargo before trains leave. – Their goal: move towards a fixed schedule. – Used OR techniques to determine the best fixed schedule – Savings in cost: $170 million. 41 The 2003 Edelman Prize for OR/MS Practice "Planning the UPS Air Network." Goal: develop an optimize service network design for express package delivery. – minimal-cost aircraft routes – fleet assignments – allocation of packages to routes – 17,000 origin-destination pairs, 160 aircraft z Savings from optimization: $275 million 42 Summary z Gauss-Jordan solving of equations and other background in linear algebra z Optimization in Radiation Delivery z Modeling in practice is an art form. It requires finding the right simplifications of reality for a given situation. z The 2003 Edelman prize competition. 43 z Animation of the Gauss-Jordan Elimination Algorithm 44 Solving a System of Equations x1 x2 x3 x4 1 2 -1 2 1 1 4 -1 2 1 -1 2 = = = 0 6 -3 To solve a system of equations, use Gauss-Jordan elimination. 45 To solve the system of equations: x1 x2 x3 x4 1 2 4 1 = 0 2 1 -1 -1 = 6 -1 1 2 2 = -3 46 Pivot on the element in row 1 column 1 x1 x2 x3 x4 1 2 4 1 = 0 2 0 -3 1 -1 -9 -1 -3 = 6 -1 0 1 3 2 6 2 3 = -3 Subtract 2 times constraint 1 from constraint 2. Add constraint 1 to constraint 3. 47 Pivot on the element in Row 2, Column 2 x1 x2 x3 x4 1 2 4 1 = 0 2 0 -3 1 -1 -9 -1 -3 = 6 -1 0 1 3 2 6 2 3 = -3 Take two minutes and carry out the “pivot” operation. 48 Pivot on the element in Row 2, Column 2 x1 x2 x3 x4 1 2 0 -2 4 -1 1 = 0 4 0 -3 1 -9 3 -3 1 = -2 6 0 3 0 -3 6 3 0 = -3 3 Divide constraint 2 by -3. Subtract multiples of constraint 2 from constraints 1 and 3. 49 Pivot on the element in Row 3, Column 3 x1 x2 x3 x4 1 0 -2 0 -1 = 4 2 0 1 3 0 1 = -2 1 0 0 -3 1 0 = -1 3 Divide constraint 3 by -3. Add multiples of constraint 3 to constraints 1 and 2. Suppose x4 = 0. What are x1, x2, x3? 50 The fundamental operation: pivoting x1 x2 x3 x4 a11 a12 a13 a14= = b1 = b2 = a21 a22 a23 a24 = a31 a32 a33 a34 = b3 Pivot on a23 51 Pivot on a23 ⎯a11 =a11 –a13(a21/a23) x1 x2 x3 x4 a1111 ⎯a a12 a013 a14= = b1 a123 = a24a/a 24 23 = = b2b/a2 23 a033 a34 a21a/a a/a 21 23 a22 22 23 a31 a32 = b3 What will be the next coefficient of b1? a32? of aij for i ≠2? 52