Exploring OMPR with HiGHS solver | R-bloggers

• Read original article here

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. Since OMPR can directly call HiGHS optimizer, we can solve the model and get solution as shown below. Solving the above problem results in an objective value of 5.75 and solution of (0.5, 2.25)

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. is the quantity to be shipped from plant to market (decision variable) Constraint (b) ensures that total supply from a plant is below capacity Constraint (c) ensures that demand for each market is met. 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) %
# 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.

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. if color is used, if state is colored with color . Objective (a) is to minimize the number of colors used Constraint (b) ensures that each state gets some color Constraint (c) ensures that if state and are adjacent, they don’t get the same color. # 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]