Exploring OMPR with HiGHS solver

There is a class of software for modeling optimization problems referred to as algebraic modeling systems which provide a unified interface to formulate optimization problems in a manner that is close to mathematical depiction and have the ability to link to different types of solvers (sparing the user from solver specific ways of formulating the problem). Both commercial and open source options are available. GAMS and AMPL are examples of commercial options. The popular open source options are JuMP in Julia and Pyomo in python. I have typically used Pyomo in Python but have explored using it from R. I recently became aware of algebraic modeling system in R provided by OMPR package developed by Dirk Schumacher.

There are commercial and open-source options available for solvers also. For a class of optimization problems referred to as Mixed Integer Linear Programs (MILP), the commercial solvers such as CPLEX, and GUROBI perform significantly better than open source solvers such as glpk, and CBC. A new open-source solver HiGHS has been developed recently that has generated quite a bit of buzz and by different accounts looks like a promising option. There is now a highs package in R that can call the HiGHS solver.

In this blog, I wanted to explore using OMPR modeling system with HiGHS solver by using it to solve a few examples of LP/MILP problems.

Example 1: Example from highs package

Here I want to just describe the example in mathematical notation and show how OMPR model is close to mathematical notation. The full details of this example are in this location.

mdl = MIPModel() %>%
      add_variable(x0, lb = 0, ub = 4, type = "continuous") %>%
      add_variable(x1, lb = 1, type = "continuous") %>%
      set_objective(x0+x1+3, sense = "min") %>%
      add_constraint(x1 <= 7) %>%
      add_constraint(x0 + 2*x1 <= 15) %>%
      add_constraint(x0 + 2*x1 >= 5) %>%
      add_constraint(3*x0 + 2*x1 >= 6)

Since OMPR can directly call HiGHS optimizer, we can solve the model and get solution as shown below.

# solve model
s = mdl %>% solve_model(highs_optimizer())

# get solution
s$status
s$objective_value
s$solution

Solving the above problem results in an objective value of 5.75 and solution of (0.5, 2.25)

Example 2: Transportation Problem

This example discusses a transporation problem from GAMS model library where the goal is to find the minimum cost way to meet market demand with available plant capacity. We just show how the OMPR package can handle variables involving indices using this example. The full description of this example is in this location.

np = length(plants)
nm = length(mkts)
# create ompr model
mdl = MIPModel() %>%
  add_variable(x[i, j], i=1:np, j=1:nm, type = "continuous",lb = 0) %>%
  # objective: min cost
  set_objective(sum_over(cost(i, j) * x[i, j], i = 1:np, j = 1:nm), sense = "min") %>% 
  # supply from each plant is below capacity
  add_constraint(sum_over(x[i, j], j = 1:nm) <= cap[i], i = 1:np) %>%  
  # supply to each market meets demand
  add_constraint(sum_over(x[i, j], i = 1:np) >= dem[j], j = 1:nm)

The figure on the left show the supply network (plants on top and markets below with numbers being capacity for plants and demand for markets). The figure on the right shows the solution where Chicago market is supplied by Seattle plant and San Diego plant supplies both New York and Topeka markets.

Example 3: Map Coloring Problem

This example discusses a map coloring problem where the goal is to use the minimum number of colors so that no two adjacent states in the US map have the same color. In this example also, I am just showing the mathematical formulation and OMPR model. The full description of this example is in this location.

# OMPR model
ns = nrow(nodes_df)
nc = 4
edge_str = edge_df %>% mutate(edge_str = glue("{fromid}_{toid}")) %>% pull(edge_str)
mdl = MIPModel()
mdl = mdl %>% add_variable(x[i, c], i = 1:ns, c = 1:nc, type = "integer", lb = 0, ub = 1)
mdl = mdl %>% add_variable(y[c], c = 1:nc, type = "integer", lb = 0, ub = 1)
mdl = mdl %>% set_objective(sum_over(y[c], c=1:nc), sense = "min")
mdl = mdl %>% add_constraint(sum_over(x[i, c], c = 1:nc) == 1, i = 1:ns)
mdl = mdl %>% add_constraint(x[i, c] + x[j, c] <= y[c], i = 1:ns, j = 1:ns, c = 1:nc, glue("{i}_{j}") %in% edge_str)

Solving this problem give the following map coloring