Optimization: simply do more with less, zoo, buses and kids (Part2, python, java, C++…)

Languages bar chart generated from OPL CPLEX through a python call

Part 1 is now in 23 human languages, 23 like the number of chromosomes we humans all have and now it’s time to give more options with computer languages.

IBM ILOG CPLEX Optimization Studio is a prescriptive analytics solution that enables rapid development and deployment of decision optimization models using mathematical and constraint programming.

In that previous post, I gave two examples with the most popular languages for CPLEX : OPL and Python, an algebraic modeling language (AML), and a general programming language (GPL)

Here, I will show many more computer options than OPL and Python docplex (which are the two most common ones).

API pie chart generated from OPL CPLEX through a Python call

With both OPL and Python docplex API, we can

• call CPLEX Linear Programming and Constraint Programming (CPOptimizer)
• call CPLEX on a given machine, in the cloud (WML) and on an optimization server

I really recommend OPL and the use of OPL C++/JAVA/.NET/Python API. For the ones who really need to write Python only, I suggest the Python docplex api.

But there are many other options with CPLEX.

1. docplex is not the only CPLEX Python API
2. JAVA, C++, Scala, Julia, Optano
3. CPLEX specific formats
4. Other AML (Algebraic Modeling Languages) like GAMS , AIMMS and AMPL

But also with other tools that won’t scale

1. opensource like GLPK
2. geometry
3. plain enumeration

In the past, this bus and kids example has not only helped many to understand what optimization is, but also some start with a small model and then improve. So let me share what helped me.

python docplex

from docplex.mp.model import Model

mdl = Model(name='buses')
nbbus40 = mdl.integer_var(name='nbBus40')
nbbus30 = mdl.integer_var(name='nbBus30')
mdl.add_constraint(nbbus40*40 + nbbus30*30 >= 300, 'kids')
mdl.minimize(nbbus40*500 + nbbus30*400)

mdl.solve(log_output=True,)

for v in mdl.iter_integer_vars():
    print(v," = ",v.solution_value)

docplex is not the only CPLEX python API

CPLEX matrix python API

# Import packages.
import cplex
my_prob = cplex.Cplex()
my_obj = [500, 400]
my_ub = [cplex.infinity, cplex.infinity]
my_lb = [0.0, 0.0]
my_ctype = "II"
my_colnames = ["nbBus40", "nbBus30"]
my_rhs = [300]
my_rownames = ["nbKids"]
my_sense = "G"
my_prob.objective.set_sense(my_prob.objective.sense.minimize)
my_prob.variables.add(obj=my_obj, lb=my_lb, ub=my_ub, types=my_ctype,
                       names=my_colnames)
rows = [[["nbBus40", "nbBus30"], [40,30]]]
my_prob.linear_constraints.add(lin_expr=rows, senses=my_sense,
                                rhs=my_rhs, names=my_rownames)
my_prob.solve()
print("cost  = ", my_prob.solution.get_objective_value())
numcols = my_prob.variables.get_num()
    
sol = my_prob.solution.get_values()
for j in range(numcols):
    print(my_colnames[j],"  = ", sol[j])

Pyomo

import pyomo.environ as pyo
from pyomo.opt import SolverFactory
opt = pyo.SolverFactory('cplex')
model = pyo.ConcreteModel()
model.nbBus = pyo.Var([40,30], domain=pyo.PositiveIntegers)
model.OBJ = pyo.Objective(expr = 500*model.nbBus[40] + 400*model.nbBus[30])
model.Constraint1 = pyo.Constraint(expr = 40*model.nbBus[40] + 30*model.nbBus[30] >= 300)
results = opt.solve(model)
print("nbBus40=",model.nbBus[40].value)
print("nbBus30=",model.nbBus[30].value)

Pulp

import pulp
import cplex
bus_problem = pulp.LpProblem("bus", pulp.LpMinimize)
nbBus40 = pulp.LpVariable('nbBus40', lowBound=0, cat='Integer')
nbBus30 = pulp.LpVariable('nbBus30', lowBound=0, cat='Integer')
# Objective function
bus_problem += 500 * nbBus40 + 400 * nbBus30, "cost"
# Constraints
bus_problem += 40 * nbBus40 + 30 * nbBus30 >= 300
bus_problem.solve(pulp.CPLEX())
print(pulp.LpStatus[bus_problem.status])
for variable in bus_problem.variables():
    print ("{} = {}".format(variable.name, variable.varValue))

CVXPY

# Import packages.
import cvxpy as cp
# Define and solve the CVXPY problem.
nbBus40 = cp.Variable(integer=True)
nbBus30 = cp.Variable( integer=True)
cost = 500*nbBus40+400*nbBus30
prob = cp.Problem(cp.Minimize(cost),[40*nbBus40+30*nbBus30>=300,
                                     nbBus40>=0,nbBus30>=0
                                     ])
prob.solve(solver=cp.CPLEX,verbose=True)
# Print result.
print("\nThe minimal cost is", prob.value)
print("number buses 40 seats = ",nbBus40.value)
print("number buses 30 seats = ",nbBus30.value)

JAVA, C++, Scala,Julia

JAVA

import ilog.concert.*;
import ilog.cplex.*;
import ilog.opl.*;
static public void main(String[] args) throws Exception
{
     IloCplex cplexBus = new IloCplex();
// decision variables
     IloNumVar nbbus40 = cplexBus.numVar(0, 10000,IloNumVarType.Int); 
     IloNumVar nbbus30 = cplexBus.numVar(0, 10000,IloNumVarType.Int);
// move at least 300 kids to the zoo
     cplexBus.addGe(cplexBus.sum(cplexBus.prod(40,nbbus40), cplexBus.prod(30,nbbus30)),300);
// objective : minimize cost = 500*nbbus40+400*nbBus30
     cplexBus.addMinimize(cplexBus.sum(cplexBus.prod(500,nbbus40), cplexBus.prod(400,nbbus30)));
cplexBus.solve();
System.out.println("nbbus40 : "  +cplexBus.getValue(nbbus40) );
      System.out.println("nbbus30 : "  +cplexBus.getValue(nbbus30) );
}

C++

#include <iostream>
#include <ilcplex/ilocplex.h>
using std::cout;
using std::cerr;
using std::endl;
int
main(void)
{
   IloEnv env;
   try {
      IloModel model(env);
// int nbKids = 300
      int nbKids = 300;
      // float costBus40 = 500
      // float costBus30 = 400
      double costBus40 = 500;
      double costBus30 = 400;
// dvar int+ nbBus40;
      // dvar int+ nbBus30;
      IloIntVar nbBus40(env, 0, IloInfinity, "nbBus40");
      IloIntVar nbBus30(env, 0, IloInfinity, "nbBus30");
// minimize costBus40 * nbBus40 + costBus30 * nbBus30
      IloExpr cost = costBus40 * nbBus40 + costBus30 * nbBus30;
      model.add(IloMinimize(env, cost));
// subjec to {
      //  40 * nbBus40 + 30 * nbBus30 >= nbKids
      // }
      model.add(40 * nbBus40 + 30 * nbBus30 >= nbKids);
IloCplex cplex(model);
      cplex.solve();
cout << "Use " << cplex.getValue(nbBus40) << " buses of type 40" << endl;
      cout << "Use " << cplex.getValue(nbBus30) << " buses of type 30" << endl;
      cout << "Total cost: " << cplex.getValue(cost) << endl;
env.end();
   }
   catch (IloException &e) {
      cerr << "exception: " << e << endl;
      throw;
   }
   return 0;
}

SCALA

import com.decisionbrain.cplexscala.Modeler._

object Buses {

  def solve(): Unit = {

    val nbKids: Int = 300
    val costBust40: Double = 500
    val costBust30: Double = 400

    implicit val model: MpModel = MpModel("buses")

    val nbBus40 = model.intVar(name="nbBus40")
    val nbBus30 = model.intVar(name="nbBus30")

    model.add(nbBus40*40 + nbBus30*30 >= nbKids)

    model.add(minimize(costBust40*nbBus40 + costBust30*nbBus30))

    model.solve()

    println(s"# bus of size 40: ${model.getValue(nbBus40)}")
    println(s"# bus of size 30: ${model.getValue(nbBus30)}")
  }

  def main(args: Array[String]): Unit = {
    solve()
  }
}

Julia

import Pkg
Pkg.add("JuMP")
Pkg.build("CPLEX")
Pkg.add("CPLEX")
using JuMP, CPLEX
model = Model(CPLEX.Optimizer)
@variable(model, nbbus40, lower_bound=0, integer=true)
@variable(model, nbbus30, lower_bound=0, integer=true)
# Objective - minimize cost
@objective(model, Min, 500*nbbus40+400*nbbus30)
# 300 kids at least
@constraint(model,40*nbbus40+30*nbbus30>=300)
# Solve
JuMP.optimize!(model)
# Display 
print("nbbus40=",(JuMP.value(nbbus40)))
print("\n");  
print("nbbus30=",(JuMP.value(nbbus30)))

Optano

// create a new dotnet project             // > dotnet new console             // Add OPTANO Modeling             // > dotnet add package 
optano.modeling              // add this code to project.cs as a              
var nbKids = 300;             var costBus40 = 500d;             
var costBus30 = 400d;              
var nbBus40 = new Variable("nbBus40", type: VariableType.Integer);             var nbBus30 = new Variable("nbBus30", type: VariableType.Integer);              var model = new Model();             
model.AddObjective(
new Objective(costBus40 * nbBus40 + costBus30 * nbBus30));             
model.AddConstraint(40 * nbBus40 + 30 * nbBus30 >= nbKids, "Kids"); 
var solver = new CplexSolver();                          
var solution = solver.Solve(model);
              
Console.WriteLine($"Objective {solution.ObjectiveValues.Single()}");             Console.WriteLine($"nbBus40: {nbBus40.Value}");             Console.WriteLine($"nbBus30: {nbBus30.Value}");

CPLEX specific formats

LP format

Minimize
 obj1: 500 nbBus40 + 400 nbBus30
Subject To
 c1: 40 nbBus40 + 30 nbBus30 >= 300
Bounds
      nbBus40 >= 0
      nbBus30 >= 0
Generals
 nbBus40  nbBus30 
End

MPS format

NAME          Zoo
ROWS
 N  obj1    
 G  c1      
COLUMNS
    MARK0000  'MARKER'                 'INTORG'
    nbBus40   obj1                          500
    nbBus40   c1                             40
    nbBus30   obj1                          400
    nbBus30   c1                             30
    MARK0001  'MARKER'                 'INTEND'
RHS
    rhs       c1                            300
BOUNDS
 LI bnd       nbBus40   0
 LI bnd       nbBus30   0
ENDATA

IBM CPLEX OPL

int nbKids=300;
float costBus40=500;
float costBus30=400;
 
dvar int+ nbBus40;
dvar int+ nbBus30;
 
minimize
 costBus40*nbBus40  +nbBus30*costBus30;
 
subject to
{
 40*nbBus40+nbBus30*30>=nbKids;
}

Other AML (Algebraic Modeling Languages) like GAMS , AIMMS and AMPL

GAMS

Variable
 cost;Integer Variables
 nbbus40, nbbus30;Equations
 obj, c1;obj.. cost =e= 500 * nbbus40 + 400 * nbbus30;
c1..  40 * nbbus40 + 30 * nbbus30 =g= 300;option optcr = 1e-3;
Model buses / all /;
Solve buses using MIP minimizing cost;

AIMMS

AMPL

param nbKids = 300;
param costBus40 = 500;
param costBus30 = 400;

var nbBus40 integer >= 0;
var nbBus30 integer >= 0;

minimize TotalBuses:
  costBus40*nbBus40 + costBus30*nbBus30;

subject to TravelNeeds:
  40*nbBus40 + 30*nbBus30 >= nbKids;

Minizinc

Other tools that won’t scale

Open source solutions like GLPK

var nbBus40, integer, >=0;
var nbBus30, integer, >=0;
s.t. ctKids:40*nbBus40+30*nbBus30>=300;
minimize obj: 500*nbBus40+400*nbBus30;
/* solve the problem */
solve;
/* and print its optimal solution */
printf("nbBus40 = ");
printf(nbBus40);
printf("\n");
printf("nbBus30 = ");
printf(nbBus30);
printf("\n");
end;

Geometry

Plain enumeration (In Python or any language)

nbKids=300
costBus40=500
costBus30=400
minCost=100000
nbBus40min=0
nbBus30min=0
for nbBus40 in range(0,nbKids // 40+1):
   for nbBus30 in range(0,nbKids // 30+1):
     if (40*nbBus40+30*nbBus30>=nbKids):
      newCost=costBus40*nbBus40+costBus30*nbBus30
      if newCost<minCost:
        minCost=newCost
        nbBus40min=nbBus40
        nbBus30min=nbBus30
        
print("Result")
print("nbBus40=",nbBus40min)
print("nbBus30=",nbBus30min)

or in Scratch for kids that do not study Python

Or even with IBM Quantum Computing (simulator or real quantum device)

PS:

My bus and kids example escaped me. Others started to use it:

1. JavaSolver
2. optimiser.com
3. picat

Picat> K=300,S=[30,40],P=[400,500],N=S.len,   B=new_list(2),B::0..100,sum([S[I]*B[I]:I in 1..N])#>=K,   Z#=sum([B[I]*P[I]:I in 1..2]),solve($[min(Z)],B),println([Z,B])

The goal of that article is to show many programming options with the same tiny example. (Zoo example)

3 basic concepts:

2 decision variables : number of buses 40 seats and 30 seats

1 objective : minimize cost

1 constraint : move 300 kids to the zoo

And many more information in Low barrier to entry to optimization