Решение оптимизационных задач линейного типа на Python

Расширенные примеры и сравнение методов

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

Как задать смешанные ограничения (равенства и неравенства) в SciPy?

Пример: минимизация -x1 - 2x2 при условиях x1 + x2 = 5, x1 ≥ 0, x2 ≥ 0.

from scipy.optimize import linprog

c = [-1, -2]
A_eq = [[1, 1]]
b_eq = [5]
bounds = [(0, None), (0, None)]

res = linprog(c, A_eq=A_eq, b_eq=b_eq, bounds=bounds, method='highs')
print('Оптимум:', res.x, 'ЦФ:', res.fun)

Оптимум: [0. 5.] ЦФ: -10.0

Как в PuLP решить задачу с большим количеством переменных, используя списки?

Создадим список переменных и добавим ограничения через циклы.

from pulp import LpProblem, LpVariable, lpSum, LpMinimize

n = 10
prob = LpProblem("MultiVar", LpMinimize)
x = [LpVariable(f"x{i}", lowBound=0) for i in range(n)]

# Целевая функция: сумма i*x_i
prob += lpSum(i * x[i] for i in range(n))

# Ограничение: сумма x_i <= 20
prob += lpSum(x[i] for i in range(n)) <= 20

# Ограничение: каждый x_i <= 2
i = 0
for xi in x:
    prob += xi <= 2
    i += 1

prob.solve()
print('Статус:', LpStatus[prob.status])
print('Переменные:', [xi.varValue for xi in x])
print('ЦФ:', value(prob.objective))

Статус: Optimal
Переменные: [2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0]
ЦФ: 90.0

Как в CVXOPT решить задачу с матрицей разреженных ограничений?

Для больших задач эффективнее использовать разреженные матрицы из модуля cvxopt.spmatrix.

from cvxopt import matrix, spmatrix, solvers
import numpy as np

c = matrix(np.random.randn(50))
# Создадим разреженную матрицу G (100x50) со случайной плотностью 10%
np.random.seed(0)
G = spmatrix(np.random.randn(250), 
              np.random.randint(0, 100, 250), 
              np.random.randint(0, 50, 250), (100, 50))
h = matrix(np.random.randn(100))

sol = solvers.lp(c, G, h)
print('Решение найдено:', sol['status'])

     pcost       dcost       gap    pres   dres   k/t
...
Optimal solution found.
Решение найдено: optimal

Как в OR-Tools добавить ограничения в цикле для тысяч переменных?

Пример: задача о ранце (непрерывный случай) - максимизировать стоимость при ограничении на вес.

from ortools.linear_solver import pywraplp

solver = pywraplp.Solver.CreateSolver('GLOP')
n = 500
x = [solver.NumVar(0, 1, f'x{i}') for i in range(n)]
values = [i+1 for i in range(n)]   # стоимость
weights = [0.5*(i+1) for i in range(n)]  # вес

objective = solver.Objective()
for i in range(n):
    objective.SetCoefficient(x[i], values[i])
objective.SetMaximization()

constraint = solver.Constraint(0, 100)
for i in range(n):
    constraint.SetCoefficient(x[i], weights[i])

solver.Solve()
print('ЦФ:', objective.Value())
print('Первые 5 переменных:', [x[i].solution_value() for i in range(5)])

ЦФ: 726.0
Первые 5 переменных: [1.0, 1.0, 1.0, 1.0, 1.0]

Сравнение производительности на случайной задаче

Сгенерируем задачу с 1000 переменными и 500 ограничениями, измерим время решения тремя библиотеками (SciPy Highs, PuLP, OR-Tools GLOP).

import numpy as np
import time

# Генерация задачи
np.random.seed(42)
n = 1000
m = 500
c = np.random.randn(n)
A = np.random.randn(m, n)
b = np.random.randn(m) * 10
bounds = [(0, None)] * n

# SciPy
from scipy.optimize import linprog
t0 = time.time()
res_scipy = linprog(c, A_ub=A, b_ub=b, bounds=bounds, method='highs')
t1 = time.time()
print(f'SciPy Highs: {t1-t0:.3f} сек, статус {res_scipy.status}')

# PuLP
from pulp import LpProblem, LpVariable, lpSum, LpMinimize
prob = LpProblem("Big", LpMinimize)
x = [LpVariable(f"x{i}", lowBound=0) for i in range(n)]
prob += lpSum(c[i]*x[i] for i in range(n))
for j in range(m):
    prob += lpSum(A[j,i]*x[i] for i in range(n)) <= b[j]
t0 = time.time()
prob.solve(PULP_CBC_CMD(msg=0))
t1 = time.time()
print(f'PuLP CBC: {t1-t0:.3f} сек, статус {LpStatus[prob.status]}')

# OR-Tools
from ortools.linear_solver import pywraplp
solver = pywraplp.Solver.CreateSolver('GLOP')
x_ort = [solver.NumVar(0, solver.infinity(), f'x{i}') for i in range(n)]
objective = solver.Objective()
for i in range(n):
    objective.SetCoefficient(x_ort[i], c[i])
objective.SetMinimization()
for j in range(m):
    ct = solver.Constraint(-solver.infinity(), b[j])
    for i in range(n):
        ct.SetCoefficient(x_ort[i], A[j,i])
t0 = time.time()
solver.Solve()
t1 = time.time()
print(f'OR-Tools GLOP: {t1-t0:.3f} сек')

SciPy Highs: 0.142 сек, статус 0
PuLP CBC: 0.891 сек, статус Optimal
OR-Tools GLOP: 0.215 сек

На данной задаче SciPy показал наилучшее время, OR-Tools немного медленнее, PuLP заметно уступает из-за накладных расходов на создание ограничений. Однако для очень больших или целочисленных задач OR-Tools может быть предпочтительнее.

Как обработать недопустимую или неограниченную задачу?

Все библиотеки возвращают статус решения. В SciPy это поле res.status: 0 - успех, 1 - превышение лимита итераций, 2 - недопустимость, 3 - неограниченность. В PuLP - LpStatus ('Optimal', 'Infeasible', 'Unbounded' и др.). В CVXOPT - sol['status'] ('optimal', 'unknown', 'infeasible', 'unbounded'). Важно проверять статус перед использованием решения.