MITx 15.071x - Unit 8 - Linear Optimization - Radiation Therapy: An Application of Linear Optimization

linear-optimization.rb

# Unit 8 of the Analytics Edge
# 
# Minimize : 
#   z = (1 + 2) * X1 + (2 + 2.5) * X2 + 2.5 * X3 + X4 + 2 * X5 + (1 + 2 + 1) * X6
#
# Beamlet 
# 1 - 1   2   2   0   0   0   0   0   0
# 2 - 0   0   0   1   2   2.5 0   0   0
# 3 - 0   0   0   0   0   0   1.5 1.5 2.5
# 4 - 1   0   0   2   0   0   1   0   0
# 5 - 0   1   0   0   2   0   0   1   0
# 6 - 0   0   1   0   0   2   0   0   1
#
# Constraints:
# p - voxel at number 2 gets a minimum dose of 7
# q -                 4
# r -                 7
# s -                 8
# t -                 5 gets a maximal dose of 5

require 'rglpk'
require 'matrix'

p = Rglpk::Problem.new
p.name = "The Analytics Edge - MITx 15.071x - Unit 8 - Linear Optimization - Radiation Therapy: An Application of Linear Optimization"
p.obj.dir = Rglpk::GLP_MIN

rows = p.add_rows(5)
rows[0].name = "p"
rows[0].set_bounds(Rglpk::GLP_LO, 7, 0)
rows[1].name = "q"
rows[1].set_bounds(Rglpk::GLP_LO, 7, 0)
rows[2].name = "r"
rows[2].set_bounds(Rglpk::GLP_LO, 7, 0)
rows[3].name = "s"
rows[3].set_bounds(Rglpk::GLP_LO, 7, 0)
rows[4].name = "t"
rows[4].set_bounds(Rglpk::GLP_UP, 0, 5)

cols = p.add_cols(6)
cols[0].name = "x1"
cols[0].set_bounds(Rglpk::GLP_LO, 0.0, 0.0)
cols[1].name = "x2"
cols[1].set_bounds(Rglpk::GLP_LO, 0.0, 0.0)
cols[2].name = "x3"
cols[2].set_bounds(Rglpk::GLP_LO, 0.0, 0.0)
cols[3].name = "x4"
cols[3].set_bounds(Rglpk::GLP_LO, 0.0, 0.0)
cols[4].name = "x5"
cols[4].set_bounds(Rglpk::GLP_LO, 0.0, 0.0)
cols[5].name = "x6"
cols[5].set_bounds(Rglpk::GLP_LO, 0.0, 0.0)

radiation_matrix = Matrix[ 
  [ 1,   2,   2,   0,   0,   0,   0,   0,   0 ], 
  [ 0,   0,   0,   1,   2,   2.5, 0,   0,   0],
  [ 0,   0,   0,   0,   0,   0,   1.5, 1.5, 2.5],
  [ 1,   0,   0,   2,   0,   0,   1,   0,   0],
  [ 0,   1,   0,   0,   2,   0,   0,   1,   0],
  [ 0,   0,   1,   0,   0,   2,   0,   0,   1]
]

healthy_matrix = Matrix[ [1], [0], [1], [0], [1], [1], [0], [0], [1] ] # 1 is healthy, 0 is tumor

coefs = (radiation_matrix * healthy_matrix).column(0).to_a
p.obj.coefs = coefs

p.set_matrix(
radiation_matrix.column(2 - 1).to_a +
radiation_matrix.column(4 - 1).to_a + 
radiation_matrix.column(7 - 1).to_a + 
radiation_matrix.column(8 - 1).to_a +
radiation_matrix.column(5 - 1).to_a
)

p.simplex
z = p.obj.get
x1 = cols[0].get_prim
x2 = cols[1].get_prim
x3 = cols[2].get_prim
x4 = cols[3].get_prim
x5 = cols[4].get_prim
x6 = cols[5].get_prim

puts "z: #{z}"
puts "x1: #{x1}"
puts "x2: #{x2}"
puts "x3: #{x3}"
puts "x4: #{x4}"
puts "x5: #{x5}"
puts "x6: #{x6}"