SobesNBI

The solution of some integer optimization problem described here.

Files

The main files are:

• crash_test.py - crash test
• LICENSE - MIT license
• main.py - main script
• problem.py - problem description
• README.md - this file
• requirements.txt - requirements file for pip install
• sample_problem.csv - sample problem
• solver.py - solver
• test_solver.py - manually written tests for solver
• Цифровой планировщик.pdf - the original problem solved in this project

The following files relate to 2nd part of task:

• crash_test_ext.py - crash test for extended problem solver
• grammar_parser.py - the parser for domain specific language describing the problem (for part 2 of the task)
• main_ext.py - main script for part 2 of the task
• problem_ext.py - extended problem description (for part 2 of the task)
• problem.schedlang - sample problem written on DSL (domain specific language)
• solver_ext.py - solver for extended problem

Install the dependencies

1. Create the virtual environment.
2. Activate it
3. Install requirements from requirements.txt

The following code do the things mentioned above:

python3 -m venv .env
source .env/bin/activate
pip install -r requirements.txt

How to run

Just run from the command line:

python3 main.py <filename.csv>

e.g.

python3 main.py sample_problem.py

This script reads <filename.csv>, solves the problem and prints the solution. If this script is run without arguments it generates and solves some random problem. The generated problem is saved to random_problem.csv.

The csv-formatted table should have the following columns:

• id is the id of the product (integer)
• t is the time required to process product,
• d is the deadline,
• p is the prize (reward) if the product is ready before deadline.

The first line is treated as the header and ignored. The number of products equals the number of rows.

Below is an example of csv-file:

 id,            time,        deadline,          reward
  1,    0.4170220047,     2.424675303,    0.1403869386
  2,    0.7203244934,     1.886375315,    0.1981014891
  3, 0.0001143748173,     3.083487752,    0.8007445687
  4,    0.3023325726,    0.9200351238,    0.9682615757
  5,    0.1467558908,     3.951528464,    0.3134241782
  6,   0.09233859477,    0.1232441694,    0.6923226157
  7,    0.1862602114,     3.017103796,    0.8763891523
  8,     0.345560727,     1.877871611,    0.8946066635
  9,    0.3967674742,     2.514104228,   0.08504421137

The solution is printed on the screen as follows:

Schedule (solution):
[6, 1, 4, 8, 2, 5, 9, 3, 7]
Reward 4.96928139317 out of 4.96928139317

The schedule shows the order in which products are processed. Reward is a sum of all prizes got for all products finished before the deadline has expired. The second number is the total unconditional sum of prizes for all products. This plays a role of the upper bound limit for the actual reward. In the case above, the actual reward is equal to the max possible reward.

How to run tests

Just run

pytest

with no arguments. These tests are generated by-hands and check solver on the set of problems with known solutions.

How to run crash test

Crash test is a sequence of tests with randomly generated problems. The total reward found by pulp-based solver is compared with reward found by the brute-force solver. The brute-force solver is based on the enumeration of all possible schedules and then choosing the best one.

The crash test can be run with the following command:

python3 crash_test.py

Extensions

The extensions corresponds to the part 2 of the original task (see Цифровой планировщик.pdf). They are represented in files with ext suffix. These file contains ExtendedProblem and ExtendedPulpSolver classes which extends Problem and PulpSolver classes with constraints of a new type and fines. The constraint of the new type requires that some product should be processed before another one. The fines are subtracted from the total reward if the product is got ready after the deadline has expired.

This module also includes parser for Backus-Naur form domain-specific language (BNF DSL) that is suited to describe problems with extended constraints and fines. The parser is represented by class GrammarParser in grammar_parser.py.

To run extended module, print

python3 main_ext.py

It will read and parse problem.schedlang file, run extended pulp-based solver and compare the solution with the brute force solution.

To run crash tests, print

python3 crash_test_ext.py

Schedlang

The example of problem is presented below

# format: prod_name = t, d, p
# prod_name - is arbitrary name for product (letters + digits, starts from letter)
# t - duration, d - deadline, p - award
product prod0 = 0.417022004702574, 2.424675303015106, 0.14038693859523377
product prod1 = 0.7203244934421581, 1.8863753148148266, 0.1981014890848788
product prod2 = 0.00011437481734488664, 3.083487751785418, 0.8007445686755367
product prod3 = 0.30233257263183977, 0.9200351237918285, 0.9682615757193975
product prod4 = 0.14675589081711304, 3.9515284637592543, 0.31342417815924284
product prod5 = 0.0923385947687978, 0.12324416939066773, 0.6923226156693141
product prod6 = 0.1862602113776709, 3.01710379580281, 0.8763891522960383
product prod7 = 0.34556072704304774, 1.8778716106520714, 0.8946066635038473
product prod8 = 0.39676747423066994, 2.5141042280058823, 0.08504421136977791

# format: prod_name_1 before prod_name_1
prod2 before prod3
prod8 before prod0

# format: fine for prod_name is f
fine for prod0 is 1.0
fine for prod2 is 1
fine for prod5 is 10

Developers notes

What variables are used

• logical variables x[i,j] encodes the schedule, if x[i,j] = 1 then the product j is placed on the i-th step in the schedule. Constrains: each row and each column of array x should have only one non-zero element.
• Rows u[:,j] encode if j-th product is finished before the deadline. If it is so, u[i,j] = 1, where i is the column number for which x[i,j] = 1.
• Rows v[:,j] encode if j-th product is finished after the deadline has expired and a fine is applied.

Variables u and v can be expressed as step-functions:

$\begin{align} u_{i,j} &= x_{i,j}*h(d_j - \tau_i), \\ v_{i,j} &= x_{i,j}*h(\tau_i - d_j) \end{align}$

where $d_j$ is a deadline for $j$-th product, $\tau_i$ is time when $i$-th step of the schedule will be completed. It is calculated as follows

${\tau_i = \sum\limits_{k=0}^{i}\sum\limits_{j=0}^{n-1} {x_{i,j}t_j,}}$

where $n$ is the number of products, $t_j$ is time required to process product $i$. The next subsection explains how to implement the step function.

How to implement Heaviside function if solver can solve only linear problem

If the objective function involves calculation of the Heaviside step function

$h(\xi) = \begin{cases} 1, \textrm{if }\xi\geq 0\\ 0, \textrm{otherwise}. \end{cases}$

it becomes non-linear. Since mixed integer programming works only with linear function, we need to transform the original non-linear problem to linear one. This can be done in the following way.

Let the objective function is defined as

$f = \sum\limits_{i}^{n} h(\xi_i) p_i$,

where $\xi_i$ are some variables and $p_i$ are some constant. We define the new variables $u_i = h(\xi_i)$ and add two constraints per each new variable:

$ -M(1-u_i) \leq \xi_i < Mu_i$

The constant $M$ should satisfy the inequality: $M > \sup|\xi_i| + \varepsilon$, where $\varepsilon$ is some small positive number. These constraints are satisfied if and only if $u_i = h(\xi_i)$ holds. Indeed, if $x \geq 0$ the right-hand side plays, if $x < 0$ the left-hand side plays.

How to choose value for M

M can have any value greater than $\sup|\xi_i|$. But if we take some value like 1e+10 the solver will produce invalid results: the logical variables x[i,j] becomes fractional. I assume that this is due to truncated precision of floating point values when the they are passed to the selected backend solver in the text format. The values are truncated to 13 digits (see mps-files written by PuLP in /tmp/ folder on your computer). This can be fixed (16 digits output will give max precision) but this problem will arise again due to the limited precision of doubles itself. So, in this project a special attention was paid to the values of M (estimation of $sup|\xi_i|$).

Why GUROBI solver is used instead of default one

It seems that default CBC solver v. 2.10.3 installed by the latest at this moment PuPL v.3.3.0 installer is broken. The CBC solvers 2.10.6 and above up to 2.10.12 are OK. To update solver manually download the binary for CBC solver and update file .env/lib/python3.12/site-packages/pulp/solverdir/cbc/linux/i64/cbc with it.

In this project, we simply switch the CBC solver to GUROBI.