Construction site

Construction site#

Click –> Live Code to activate live coding on this page! Unfortunately, the pymoo package doesn’t work yet in this book. Download this page to run the cells on your own computer.

Problem#

../_images/construction_site.jpeg

Fig. 5.1 Created using DALL-E 3#

As managers of a construction site, we have to decide the required equipment to move a given volume of excavated soil. In order to achieve the deadline for this working unit, we need to guarantee an averaged efficiency of \(2700\) \(\cfrac{\text{m}^3}{h}\) during one month.

We have our own equipment, and also we can sub-contract another (just one) company. The efficiency of each equipment and the cost are given in the table:

\( \def\euro{\unicode{x20AC}} \)

Own equip.

Company 1

Company 2

Eff. \(\left(\frac{\text{m}^3}{h}\right)\)

Cost \(\left(\frac{\euro}{h}\right)\)

Eff. \(\left(\frac{\text{m}^3}{h}\right)\)

Cost \(\left(\frac{\euro}{h}\right)\)

Eff. \(\left(\frac{\text{m}^3}{h}\right)\)

Cost \(\left(\frac{\euro}{h}\right)\)

200

500

470

4.000

640

5.400

240

800

700

5.700

730

5.500

265

1.000

800

6.500

775

6.800

330

1.500

Then, we can use part or all of our own equipment with some of the equipment options provided by another company.

What is the optimal equipment combination that minimizes the cost?

Model#

We’ll define the model as follows:

  • Design variables: Whether or not to use each type of equipment

  • Objective function: sum of the costs of the selected equipment

  • Constraint functions: make use of one sub-contractor and create an efficiency of at least \(2700\) \(\cfrac{\text{m}^3}{h}\)

Design variables#

Let’s start with our design variables. let’s define them as a list of integer values either being \(0\) or \(1\), which selects the equipment type.

(5.1)#\[\begin{split}x=\left[ \begin{matrix} {{x}_{\text{First item of own equipment}}} \\ {{x}_{\text{Second item of own equipment}}} \\ {{x}_{\text{Third item of own equipment}}} \\ {{x}_{\text{Fourth item of own equipment}}} \\ {{x}_{\text{First item of company 1}}} \\ {{x}_{\text{Second item of company 1}}} \\ {{x}_{\text{Third item of company 1}}} \\ {{x}_{\text{First item of company 2}}} \\ {{x}_{\text{Second item of company 2}}} \\ {{x}_{\text{Third item of company 2}}} \\ \end{matrix} \right]=\left[ \begin{matrix} {{x}_{1}} \\ {{x}_{2}} \\ {{x}_{3}} \\ {{x}_{4}} \\ {{x}_{5}} \\ {{x}_{6}} \\ {{x}_{7}} \\ {{x}_{8}} \\ {{x}_{9}} \\ {{x}_{10}} \\ \end{matrix} \right]\end{split}\]

The design variables either being \(0\) or \(1\) can be defined mathemetically:

(5.2)#\[x_i\in \left\{ 0,1 \right\}\text{ } i=1,2,...,10\]

Objective function#

Now we can define the objective function as the product of the dimension to represent \(\mathop {\min }\limits_x f\left(x\right) \) in (2.1):

(5.3)#\[f\left(x\right) = 500 \cdot x_1 + 800 \cdot x_2 + 1000 \cdot x_3 + 1500 \cdot x_4 + 4000 \cdot x_5 + 5700 \cdot x_6 + 6500 \cdot x_7 + 5400 \cdot x_8 + 5500 \cdot x_9 + 6800 \cdot x_{10}\]

As this is a linear relation, this can be converted to matrix notation. In case of the design variables defined as in (2.2) in the form \(\mathop {\min }\limits_x {c^T}x\) with \(c\):

(5.4)#\[\begin{split}c = {{\left[ \begin{matrix} 500 & 800 & 1000 & 1500 & 4000 & 5700 & 6500 & 5400 & 5500 & 6800 \\ \end{matrix} \right]}^{T}}\end{split}\]

Constraint function#

We need to define two constraint functions, let’s start with the inequality constraint functions as \({{g}}\left(x\right) \le 0\) as defined in (2.1):

(5.5)#\[{{g}}\left(x\right) = -200 \cdot x_1 - 240 \cdot x_2 - 265 \cdot x_3 - 330 \cdot x_4 - 470 \cdot x_5 - 700 \cdot x_6 - 800 \cdot x_7 - 640 \cdot x_8 - 730 \cdot x_9 - 775 \cdot {x_{10}} +2700\]

As this is a linear relation, this can be converted to matrix notation as defined in (2.2) in the form \({A_{ub}}x \le {b_{ub}}\):

(5.6)#\[\left[ \begin{matrix} -200 & -240 & -265 & -330 & -470 & -700 & -800 & -640 & -730 & -775 \end{matrix} \right] x \le -2700 \]

Let’s add the equality constraint function. As soon as we hire one piece of equipment this should be violated. It is defined as \({{h}}\left(x\right) \le 0\) as in (3.1):

(5.7)#\[{{h}}\left(x\right) = \max \left( x_5, x_6, x_7 \right) \cdot \max \left( x_8, x_9, x_{10} \right)\]

Clearly, this function is not linear, so we cannot define it in matrix formulation.

Method#

Now let’s solve this problem using both scipy and pymoo.

Import libraries#

import numpy as np

SciPy#

For SciPy this is easy:

import scipy as sp

pymoo#

For pymoo we need to import a lot more and be more specific:

from pymoo.problems.functional import FunctionalProblem as FunctionalProblem
from pymoo.algorithms.soo.nonconvex.ga import GA
from pymoo.optimize import minimize
from pymoo.operators.sampling.rnd import IntegerRandomSampling 
from pymoo.operators.crossover.sbx import SBX
from pymoo.operators.mutation.pm import PM
from pymoo.operators.repair.rounding import RoundingRepair

import scipy as sp 
import numpy as np

Define variables#

This is done differently in SciPy and pymoo.

SciPy#

As before, we don’t need to specify our variable \(x\) itself as defined in (5.1). However, we do need to specify that we have integers (specifically only \(0\) and \(1\) which will be covered by the bounds) In SciPy we use an array of booleans to specify which design variables are integers:

integers=np.array([True,True,True,True,True,True,True,True,True,True])

pymoo#

In pymoo we don’t have to specify that have integers, but we have to adapt all of our functions later on so that the outputs are integers. Furthermore, we must specify explicitly how many design variables we have:

n_var = 10

Define bounds#

Now let’s continue with the bounds, specified by (5.1) too:

SciPy#

The bounds are defined as before:

bounds = [[0,1],
          [0,1],
          [0,1],
          [0,1],
          [0,1],
          [0,1],
          [0,1],
          [0,1],
          [0,1],
          [0,1]]

pymoo#

In pymoo, the bounds are defined as separate arrays:

xl = np.array([0,0,0,0,0,0,0,0,0,0])
xu = np.array([1,1,1,1,1,1,1,1,1,1])

Define objective function#

Let’s define the objective function as defined in (5.4).

def obj(x):
    return np.array([500, 800, 1000, 1500, 4000, 5700, 6500, 5400, 5500, 6800])@x

Define constraint function#

Let’s define the constraint function as defined in (5.6) and (5.7)

g = np.array([-200, -240, -265, -330, -470, -700, -800, -640, -730, -775])

def nonlinconfun(x):
    return max(x[4:7]) * max(x[7:])

SciPy#

In SciPy we need to store the functions in a scipy-object including the bounds

lincon = sp.optimize.LinearConstraint(g, lb=-np.inf, ub=-2700)
nonlincon = sp.optimize.NonlinearConstraint(nonlinconfun, 0, 0)

pymoo#

In pymoo the functions can be inserted directly in the problem object later on. However, this requires a function for the linear constraints instead of just the array

def linconfun(x):
    return g@x + 2700

Solve the problem#

Now let’s solve the problem using both SciPy and pymoo.

SciPy#

result_scipy = sp.optimize.differential_evolution(func = obj,bounds = bounds,constraints=[lincon,nonlincon], integrality = integers)
print(result_scipy)

             message: Optimization terminated successfully.
             success: True
                 fun: 19000.0
                   x: [ 1.000e+00  1.000e+00  0.000e+00  1.000e+00
                        1.000e+00  1.000e+00  1.000e+00  0.000e+00
                        0.000e+00  0.000e+00]
                 nit: 34
                nfev: 1226
          population: [[ 1.000e+00  1.000e+00 ...  0.000e+00  0.000e+00]
                       [ 1.000e+00  1.000e+00 ...  0.000e+00  0.000e+00]
                       ...
                       [ 1.000e+00  1.000e+00 ...  0.000e+00  0.000e+00]
                       [ 1.000e+00  1.000e+00 ...  0.000e+00  0.000e+00]]
 population_energies: [ 1.900e+04  1.900e+04 ...  1.900e+04  1.900e+04]
              constr: [array([ 0.000e+00]), array([ 0.000e+00])]
    constr_violation: 0.0
               maxcv: 0.0

Pymoo#

In pymoo we have to define the problem and algorithm as objects, and call them from the pymoo.minimize function. For the GA functions, we repair some stuff so that we only deal with integer values

problem = FunctionalProblem(n_var=n_var, objs=obj, constr_ieq=[linconfun],constr_eq = [nonlinconfun],xl=xl, xu = xu)
algorithm = GA(sampling=IntegerRandomSampling(),crossover=SBX(repair=RoundingRepair()),mutation=PM(repair=RoundingRepair()))

Now we can solve the problem

result_pymoo = minimize(problem, algorithm)
print(result_pymoo.X)
print(result_pymoo.F)

[0 1 1 1 1 1 1 0 0 0]
[19500.]
Test yourself

Questions, discussions and comments#