22.06.2020       Выпуск 340 (22.06.2020 - 28.06.2020)       Статьи

#### Линейное программирование: оптимизации

Linear programming

### Экспериментальная функция:

Ниже вы видите текст статьи по ссылке. По нему можно быстро понять ссылка достойна прочтения или нет

Просим обратить внимание, что текст по ссылке и здесь может не совпадать.

PuLP has a more convenient linear programming API than SciPy. You don’t have to mathematically modify your problem or use vectors and matrices. Everything is cleaner and less prone to errors.

Now that you have PuLP imported, you can solve your problems.

#### Example 1

You’ll now solve this system with PuLP:

The first step is to initialize an instance of `LpProblem` to represent your model:

``````# Create the model
model = LpProblem(name="small-problem", sense=LpMaximize)
``````

You use the `sense` parameter to choose whether to perform minimization (`LpMinimize` or `1`, which is the default) or maximization (`LpMaximize` or `-1`). This choice will affect the result of your problem.

Once that you have the model, you can define the decision variables as instances of the `LpVariable` class:

``````# Initialize the decision variables
x = LpVariable(name="x", lowBound=0)
y = LpVariable(name="y", lowBound=0)
``````

You need to provide a lower bound with `lowBound=0` because the default value is negative infinity. The parameter `upBound` defines the upper bound, but you can omit it here because it defaults to positive infinity.

The optional parameter `cat` defines the category of a decision variable. If you’re working with continuous variables, then you can use the default value `"Continuous"`.

You can use the variables `x` and `y` to create other PuLP objects that represent linear expressions and constraints:

>>>
``````>>> expression = 2 * x + 4 * y
>>> type(expression)
<class 'pulp.pulp.LpAffineExpression'>

>>> constraint = 2 * x + 4 * y >= 8
>>> type(constraint)
<class 'pulp.pulp.LpConstraint'>
``````

When you multiply a decision variable with a scalar or build a linear combination of multiple decision variables, you get an instance of `pulp.LpAffineExpression` that represents a linear expression.

Similarly, you can combine linear expressions, variables, and scalars with the operators `==`, `<=`, or `>=` to get instances of pulp.LpConstraint that represent the linear constraints of your model.

Note: It’s also possible to build constraints with the rich comparison methods `.__eq__()`, `.__le__()`, and `.__ge__()` that define the behavior of the operators `==`, `<=`, and `>=`.

Having this in mind, the next step is to create the constraints and objective function as well as to assign them to your model. You don’t need to create lists or matrices. Just write Python expressions and use the `+=` operator to append them to the model:

``````# Add the constraints to the model
model += (2 * x + y <= 20, "red_constraint")
model += (4 * x - 5 * y >= -10, "blue_constraint")
model += (-x + 2 * y >= -2, "yellow_constraint")
model += (-x + 5 * y == 15, "green_constraint")
``````

In the above code, you define tuples that hold the constraints and their names. `LpProblem` allows you to add constraints to a model by specifying them as tuples. The first element is a `LpConstraint` instance. The second element is a human-readable name for that constraint.

Setting the objective function is very similar:

``````# Add the objective function to the model
obj_func = x + 2 * y
model += obj_func
``````

Alternatively, you can use a shorter notation:

``````# Add the objective function to the model
model += x + 2 * y
``````

Now you have the objective function added and the model defined.

Note: You can append a constraint or objective to the model with the operator `+=` because its class, `LpProblem`, implements the special method `.__iadd__()`, which is used to specify the behavior of `+=`.

For larger problems, it’s often more convenient to use `lpSum()` with a list or other sequence than to repeat the `+` operator. For example, you could add the objective function to the model with this statement:

``````# Add the objective function to the model
model += lpSum([x, 2 * y])
``````

It produces the same result as the previous statement.

You can now see the full definition of this model:

>>>
``````>>> model
small-problem:
MAXIMIZE
1*x + 2*y + 0
SUBJECT TO
red_constraint: 2 x + y <= 20

blue_constraint: 4 x - 5 y >= -10

yellow_constraint: - x + 2 y >= -2

green_constraint: - x + 5 y = 15

VARIABLES
x Continuous
y Continuous
``````

The string representation of the model contains all relevant data: the variables, constraints, objective, and their names.

Finally, you’re ready to solve the problem. You can do that by calling `.solve()` on your model object. If you want to use the default solver (CBC), then you don’t need to pass any arguments:

`.solve()` calls the underlying solver, modifies the `model` object, and returns the integer status of the solution, which will be `1` if the optimum is found. For the rest of the status codes, see `LpStatus[]`.

You can get the optimization results as the attributes of `model`. The function `value()` and the corresponding method `.value()` return the actual values of the attributes:

>>>
``````>>> print(f"status: {model.status}, {LpStatus[model.status]}")
status: 1, Optimal

>>> print(f"objective: {model.objective.value()}")
objective: 16.8181817

>>> for var in model.variables():
...     print(f"{var.name}: {var.value()}")
...
x: 7.7272727
y: 4.5454545

>>> for name, constraint in model.constraints.items():
...     print(f"{name}: {constraint.value()}")
...
red_constraint: -9.99999993922529e-08
blue_constraint: 18.181818300000003
yellow_constraint: 3.3636362999999996
green_constraint: -2.0000000233721948e-07)
``````

`model.objective` holds the value of the objective function, `model.constraints` contains the values of the slack variables, and the objects `x` and `y` have the optimal values of the decision variables. `model.variables()` returns a list with the decision variables:

>>>
``````>>> model.variables()
[x, y]
>>> model.variables() is x
True
>>> model.variables() is y
True
``````

As you can see, this list contains the exact objects that are created with the constructor of `LpVariable`.

The results are approximately the same as the ones you got with SciPy.

Note: Be careful with the method `.solve()`—it changes the state of the objects `x` and `y`!

You can see which solver was used by calling `.solver`:

>>>
``````>>> model.solver
<pulp.apis.coin_api.PULP_CBC_CMD object at 0x7f60aea19e50>
``````

The output informs you that the solver is CBC. You didn’t specify a solver, so PuLP called the default one.

If you want to run a different solver, then you can specify it as an argument of `.solve()`. For example, if you want to use GLPK and already have it installed, then you can use `solver=GLPK(msg=False)` in the last line. Keep in mind that you’ll also need to import it:

Now that you have GLPK imported, you can use it inside `.solve()`:

``````# Create the model
model = LpProblem(name="small-problem", sense=LpMaximize)

# Initialize the decision variables
x = LpVariable(name="x", lowBound=0)
y = LpVariable(name="y", lowBound=0)

# Add the constraints to the model
model += (2 * x + y <= 20, "red_constraint")
model += (4 * x - 5 * y >= -10, "blue_constraint")
model += (-x + 2 * y >= -2, "yellow_constraint")
model += (-x + 5 * y == 15, "green_constraint")

# Add the objective function to the model
model += lpSum([x, 2 * y])

# Solve the problem
status = model.solve(solver=GLPK(msg=False))
``````

The `msg` parameter is used to display information from the solver. `msg=False` disables showing this information. If you want to include the information, then just omit `msg` or set `msg=True`.

Your model is defined and solved, so you can inspect the results the same way you did in the previous case:

>>>
``````>>> print(f"status: {model.status}, {LpStatus[model.status]}")
status: 1, Optimal

>>> print(f"objective: {model.objective.value()}")
objective: 16.81817

>>> for var in model.variables():
...     print(f"{var.name}: {var.value()}")
...
x: 7.72727
y: 4.54545

>>> for name, constraint in model.constraints.items():
...     print(f"{name}: {constraint.value()}")
...
red_constraint: -1.0000000000509601e-05
blue_constraint: 18.181830000000005
yellow_constraint: 3.3636299999999997
green_constraint: -2.000000000279556e-05
``````

You got practically the same result with GLPK as you did with SciPy and CBC.

Let’s peek and see which solver was used this time:

>>>
``````>>> model.solver
<pulp.apis.glpk_api.GLPK_CMD object at 0x7f60aeb04d50>
``````

As you defined above with the highlighted statement `model.solve(solver=GLPK(msg=False))`, the solver is GLPK.

You can also use PuLP to solve mixed-integer linear programming problems. To define an integer or binary variable, just pass `cat="Integer"` or `cat="Binary"` to `LpVariable`. Everything else remains the same:

``````# Create the model
model = LpProblem(name="small-problem", sense=LpMaximize)

# Initialize the decision variables: x is integer, y is continuous
x = LpVariable(name="x", lowBound=0, cat="Integer")
y = LpVariable(name="y", lowBound=0)

# Add the constraints to the model
model += (2 * x + y <= 20, "red_constraint")
model += (4 * x - 5 * y >= -10, "blue_constraint")
model += (-x + 2 * y >= -2, "yellow_constraint")
model += (-x + 5 * y == 15, "green_constraint")

# Add the objective function to the model
model += lpSum([x, 2 * y])

# Solve the problem
status = model.solve()
``````

In this example, you have one integer variable and get different results from before:

>>>
``````>>> print(f"status: {model.status}, {LpStatus[model.status]}")
status: 1, Optimal

>>> print(f"objective: {model.objective.value()}")
objective: 15.8

>>> for var in model.variables():
...     print(f"{var.name}: {var.value()}")
...
x: 7.0
y: 4.4

>>> for name, constraint in model.constraints.items():
...     print(f"{name}: {constraint.value()}")
...
red_constraint: -1.5999999999999996
blue_constraint: 16.0
yellow_constraint: 3.8000000000000007
green_constraint: 0.0)

>>> model.solver
<pulp.apis.coin_api.PULP_CBC_CMD at 0x7f0f005c6210>
``````

Now `x` is an integer, as specified in the model. (Technically it holds a float value with zero after the decimal point.) This fact changes the whole solution. Let’s show this on the graph:

As you can see, the optimal solution is the rightmost green point on the gray background. This is the feasible solution with the largest values of both `x` and `y`, giving it the maximal objective function value.

GLPK is capable of solving such problems as well.

#### Example 2

Now you can use PuLP to solve the resource allocation problem from above:

The approach for defining and solving the problem is the same as in the previous example:

``````# Define the model
model = LpProblem(name="resource-allocation", sense=LpMaximize)

# Define the decision variables
x = {i: LpVariable(name=f"x{i}", lowBound=0) for i in range(1, 5)}

model += (lpSum(x.values()) <= 50, "manpower")
model += (3 * x + 2 * x + x <= 100, "material_a")
model += (x + 2 * x + 3 * x <= 90, "material_b")

# Set the objective
model += 20 * x + 12 * x + 40 * x + 25 * x

# Solve the optimization problem
status = model.solve()

# Get the results
print(f"status: {model.status}, {LpStatus[model.status]}")
print(f"objective: {model.objective.value()}")

for var in x.values():
print(f"{var.name}: {var.value()}")

for name, constraint in model.constraints.items():
print(f"{name}: {constraint.value()}")
``````

In this case, you use the dictionary `x` to store all decision variables. This approach is convenient because dictionaries can store the names or indices of decision variables as keys and the corresponding `LpVariable` objects as values. Lists or tuples of `LpVariable` instances can be useful as well.

The code above produces the following result:

``````status: 1, Optimal
objective: 1900.0
x1: 5.0
x2: 0.0
x3: 45.0
x4: 0.0
manpower: 0.0
material_a: -40.0
material_b: 0.0
``````

As you can see, the solution is consistent with the one obtained using SciPy. The most profitable solution is to produce `5.0` units of the first product and `45.0` units of the third product per day.

Let’s make this problem more complicated and interesting. Say the factory can’t produce the first and third products in parallel due to a machinery issue. What’s the most profitable solution in this case?

Now you have another logical constraint: if x₁ is positive, then x₃ must be zero and vice versa. This is where binary decision variables are very useful. You’ll use two binary decision variables, y₁ and y₃, that’ll denote if the first or third products are generated at all:

`````` 1 model = LpProblem(name="resource-allocation", sense=LpMaximize)
2
3 # Define the decision variables
4 x = {i: LpVariable(name=f"x{i}", lowBound=0) for i in range(1, 5)}
5 y = {i: LpVariable(name=f"y{i}", cat="Binary") for i in (1, 3)}
6
8 model += (lpSum(x.values()) <= 50, "manpower")
9 model += (3 * x + 2 * x + x <= 100, "material_a")
10 model += (x + 2 * x + 3 * x <= 90, "material_b")
11
12 M = 100
13 model += (x <= y * M, "x1_constraint")
14 model += (x <= y * M, "x3_constraint")
15 model += (y + y <= 1, "y_constraint")
16
17 # Set objective
18 model += 20 * x + 12 * x + 40 * x + 25 * x
19
20 # Solve the optimization problem
21 status = model.solve()
22
23 print(f"status: {model.status}, {LpStatus[model.status]}")
24 print(f"objective: {model.objective.value()}")
25
26 for var in model.variables():
27     print(f"{var.name}: {var.value()}")
28
29 for name, constraint in model.constraints.items():
30     print(f"{name}: {constraint.value()}")
``````

The code is very similar to the previous example except for the highlighted lines. Here are the differences:

• Line 5 defines the binary decision variables `y` and `y` held in the dictionary `y`.

• Line 12 defines an arbitrarily large number `M`. The value `100` is large enough in this case because you can’t have more than `100` units per day.

• Line 13 says that if `y` is zero, then `x` must be zero, else it can be any non-negative number.

• Line 14 says that if `y` is zero, then `x` must be zero, else it can be any non-negative number.

• Line 15 says that either `y` or `y` is zero (or both are), so either `x` or `x` must be zero as well.

Here’s the solution:

``````status: 1, Optimal
objective: 1800.0
x1: 0.0
x2: 0.0
x3: 45.0
x4: 0.0
y1: 0.0
y3: 1.0
manpower: -5.0
material_a: -55.0
material_b: 0.0
x1_constraint: 0.0
x3_constraint: -55.0
y_constraint: 0.0
``````

It turns out that the optimal approach is to exclude the first product and to produce only the third one.