Перейти к содержанию

Python Matrix API

Начиная с релиза ArhiCloud 1.7 в рамках Python API пользователям предоставлен расширенный набор функций для создания моделей и решения задач оптимизации в матричном виде - Python Matrix API.

Python Matrix API чаще всего используется в следующих случаях:

  • Для упрощения и ускорения процесса построения моделей оптимизации, особенно в случаях, когда работа ведётся с большими объёмами данных.
  • Позволяет эффективно использовать матричные и векторные выражения, что делает его удобным инструментом для задач, требующих обработки данных в матричной форме.

Преимущества

  • Ускорение процесса моделирования. Python Matrix API значительно ускоряет процесс проектирования модели для больших данных, которые уже представлены в бинарном виде или матричной форме.
  • Экономия памяти. Эффективная работа с разреженными матрицами позволяет минимизировать использование оперативной памяти в процессе моделирования.
  • Интеграция с популярными библиотеками. MatrixAPI легко интегрируется с такими библиотеками, как NumPy и SciPy, что упрощает работу с внешними данными и расширяет функциональность.
  • Оптимизация читаемости кода. Использование матричных представлений, когда исходные данные представлены в бинарном виде или матричной форме, делает код более компактным и читаемым, а также облегчает его поддержку и модификацию.

Порядок использования

Использование Python Matrix API осуществляется аналогично основному Python API.

1. Создание объекта модели

Создание объекта модели осуществляется с использованием класса MatrixModel.

model = MatrixModel()

Модель после создания не содержит переменных или ограничений. Установка параметров, их значений осуществляется так же, как и в базовом Python API.

2. Создание переменных

Создание переменных модели осуществляется с использованием метода add_variables(lower_bounds=0.0, upper_bounds=None, obj_values=0.0, var_types, count). С помощью данного метода имеется возможность пакетного создания переменных. В качестве параметра count указывается количество создающихся переменных.

В качестве параметров lower_bounds, upper_bounds, obj_values, var_types могут быть переданы как единичные значения параметров, в том случае, если осуществляется создание нескольких единообразных переменных, так и списки значений, если значения для каждой из переменных могут различаться. Например:

add_variables(
    lower_bounds=0.0, 
    upper_bounds=300.0, 
    obj_values=1.0, 
    var_types=ap.variable_type.continuous, 
    5
)

В представленном примере осуществляется создание пяти одинаковых по характеристикам переменных. Или:

vars_type = [ap.variable_type.continuous, ap.variable_type.binary, ap.variable_type.integer]
lower_bounds = 0
upper_bounds = [100.8, ap.kArhiplexInf, 1]
obj_values = [3, 5, 2]
count = 3

model.add_variables(
    count=count, 
    lower_bounds=lower_bounds, 
    upper_bounds=upper_bounds, 
    obj_values=obj_values, 
    var_types=vars_type
)

3. Ограничения

Добавление ограничений в объект модели может быть осуществлено как отдельно каждого ограничения, так и целиком для всей задачи в матричном виде.

3.1. Добавление всех ограничений задачи с использованием матриц и векторов

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

add_constraints(A, senses, b)
  • A – двумерная матрица линейных ограничений, где количество строк соответствует числу вводимых ограничений, а количество столбцов соответствует числу переменных. Данная матрица может быть задана в формате csr (sparse row matrix) или numpy.ndarray.

  • senses – вектор типов ограничений. Этот параметр может принимать значения в виде одиночного constraint_sense, списка List[constraint_sense] или массива numpy.ndarray. Если передано только одно значение, то оно будет использоваться для всех ограничений.

  • b – вектор констант, представляющий правые части ограничений. Его можно указать в виде типа float, списка List[float] или массива numpy.ndarray. В случае, если параметр задан одним числом, оно будет использоваться для всех правых частей ограничений.

3.2. Добавление одного квадратичного или линейного ограничения

Каждое из ограничений в задачах с квадратичными ограничениями добавляется в модель отдельным вызовом метода add_constraint().

add_constraint(P=None, q=None, sense, r)
  • P – квадратная симметричная матрица, задающая квадратичную часть ограничения с равным количеством строк и столбцов, соответствующим количеству переменных. Матрица может быть представлена в формате csr (sparse row matrix) или numpy.ndarray. При считывании используется формула \(\frac{1}{2}\overline{x}^TP\overline{x}\), что обусловливает необходимость удвоения элементов главной диагонали.

  • q – вектор коэффициентов переменных, может быть задан в формате List[float], numpy.ndarray.

  • sense – тип ограничения, constraint_sense.

  • r – константное значение правой части ограничения, float.

4. Целевая функция

В зависимости от наличия или отсутствия квадратичной или линейной составляющей в целевой функции её создание может настраиваться в процессе вызова метода set_objective(P=None, q=None, sense).

  • P – квадратная симметричная матрица, задающая квадратичную часть целевой функции с равным количеством строк и столбцов, соответствующим количеству переменных. Матрица может быть представлена в формате csr (sparse row matrix) или numpy.ndarray.

  • q – вектор коэффициентов переменных, может быть задан в формате List[float], numpy.ndarray.

  • sense – тип оптимизации, objective_sense.

В случае отсутствия квадратичной части параметр P не указывается, остаётся только линейная часть, задаваемая вектором q.

5. Запуск решения задачи

Для решения задачи оптимизации после создания модели, переменных, ограничений и установки целевого критерия используется тот же метод, что и в базовом Python API - для ArhiCloud это model.solve_remote(), а для ArhiPlex - model.solve(). В качестве параметра мэппинга в данном случае передавать название файла не требуется, так как отсутствуют наименования переменных и ограничений в матричном виде.

solve_result = model.solve_remote()

В качестве результата метод возвращает результат в виде объекта класса MatrixSolveResult.

Сохранение модели в файл

Для сохранения моделей в файл также методы MatrixModel.write_lp(file_name), MatrixModel.write_mps(file_name).

Пример 1

Постановка задачи

Переменные

  • \(x_1 \in \mathbb{R}\) - непрерывная переменная, \(0 \leqslant x_1 \leqslant 100.8\),

  • \(x_2 \in \mathbb{Z}\) - целочисленная переменная, \(x_2 \geqslant 0\),

  • \(x_3 \in \mathbb{B}\) - бинарная переменная, \(x_3 \in \{0, 1\}\).

Целевая функция

\[ 3x_1 + 5x_2 + 2x_3 \to \max\limits_{x_1, x_2, x_3} \]

Ограничения

\[ 2x_1 + 4x_2 + 3x_3 \leqslant 10 \]
\[ x_1 + 3x_2 + 2x_3 \leqslant 8 \]
\[ x_1 + x_2 \leqslant 5 \]
\[ x_2 + x_3 \leqslant 3 \]
\[ x_1 + x_3 \geqslant 1 \]

Решение

import arhiplexpy as ap
import numpy as np

def get_Data():
    """
    :return: Возвращает данные для создания модели
    """

    vars_type = [ap.variable_type.continuous, ap.variable_type.binary, ap.variable_type.integer]
    lower_bounds = 0
    upper_bounds = [100.8, ap.kArhiplexInf, 1]
    obj_values = [3, 5, 2]
    count = 3

    A = np.array([
        [2, 4, 3]
        , [1, 3, 2]
        , [1, 1, 0]
        , [0, 1, 1]
        , [1, 0, 1]
    ])

    b = [10, 8, 5, 3, 1]
    senses = [
        ap.constraint_sense.less_equal
        , ap.constraint_sense.less_equal
        , ap.constraint_sense.less_equal
        , ap.constraint_sense.less_equal
        , ap.constraint_sense.greater_equal]

    return vars_type, lower_bounds, upper_bounds, obj_values, count, A, b, senses


if __name__ == "__main__":

    # Получение данных для создания модели MILP
    vars_type, lower_bounds, upper_bounds, obj_values, count, A, b, senses = get_Data()

    # Инициализация модели
    model = ap.MatrixModel()
    # Добавление переменных
    model.add_variables(
        count=count, lower_bounds=lower_bounds, upper_bounds=upper_bounds, obj_values=obj_values, var_types=vars_type
    )
    # Добавление целевой функции
    model.set_objective(sense=ap.objective_sense.maximize)
    # Добавление ограничений
    model.add_constraints(A=A, senses=senses, b=b)

    # Сохранение модели в формате LP
    model.write_lp('MILP.lp')
    # Сохранение модели в формате MPS
    model.write_mps('MILP.mps')

    # Установление параметров расчета
    model.set_dbl_param('time_limit', 10)
    model.set_bool_param('log_to_console', False)

    # Запуск оптимизации
    result = model.solve()

    # Вывод получившегося результата
    if result.get_solve_result() == ap.solve_result.success:
        if (result.get_solution_status() == ap.solution_status.optimal or \
            result.get_solution_status() == ap.solution_status.feasible):
            print(f'Status problem: {result.get_solve_result()}')
            print(f'Objective value: {result.get_objective_value()}')
            print(f'GAP: {result.get_relative_gap()}')
            print(f'TIME: {result.get_solve_time()} s')
            print(f'Solution vector: {result.get_solution_vector()}')

            # Сохранение найденного решения
            result.write_solution('result.sol')
        else:
            print(f'Status problem: {result.get_solution_status()}')
    else:
        print(f'Status process: {result.get_solve_result()}')

Пример 2

Постановка задачи

Переменные

  • \(x_1 \in \mathbb{B}\) - бинарная переменная,

  • \(x_2 \in \mathbb{Z}\) - целочисленная переменная,

  • \(x_3, x_4 \in \mathbb{R}\) - непрерывные переменные.

Целевая функция

\[ x_1^2 + 2x_2^2 + 4x_3^2 + x_4^2 + x_1 x_2 + x_3 x_4 + 3x_1 + 5.25 \to \min\limits_{x_1, x_2, x_3, x_4} \]

Ограничения

\[ 2x_1 + 3x_2 + 4x_3 + 5x_4 \leqslant 20 \]
\[ x_1 - x_2 + 2x_3 - 3x_4 \geqslant 1 \]
\[ x_1 + 2x_2 + 3x_3 + 4x_4 = 15 \]
\[ x_1 + x_2 + x_3 + x_4 \geqslant 5 \]
\[ x_1 + 3x_2 + 2x_3 + x_4 \leqslant 12 \]

Решение

import arhiplexpy as ap
import numpy as np

def get_Data():
    """
    :return: Возвращает данные для создания модели
    """

    vars_type = np.array(
        [ap.variable_type.binary, ap.variable_type.integer, ap.variable_type.continuous, ap.variable_type.continuous]
        , np.int32
    )
    lower_bounds = np.array([0, -ap.kArhiplexInf, -ap.kArhiplexInf, -ap.kArhiplexInf])
    upper_bounds = [1, ap.kArhiplexInf, ap.kArhiplexInf, ap.kArhiplexInf]
    count = 4

    # Диагональ должна быть удвоена относительно исходной матрицы!
    P = np.array([
        [2, 1, 0, 0]
        , [1, 4, 0, 0]
        , [0, 0, 8, 1]
        , [0, 0, 1, 2]
    ])
    q = np.array([3, 0, 0, 0])
    offset = 5.25

    A = np.array([
        [2, 3, 4, 5]
        , [1, -1, 2, -3]
        , [1, 2, 3, 4]
        , [1, 1, 1, 1]
        , [1, 3, 2, 1]
    ])
    b = np.array([20, 1, 15, 5, 12])
    senses = np.array(
        [
            ap.constraint_sense.less_equal
            , ap.constraint_sense.greater_equal
            , ap.constraint_sense.equal
            , ap.constraint_sense.greater_equal
            , ap.constraint_sense.less_equal
        ]
        , np.int32
    )

    return vars_type, lower_bounds, upper_bounds, count, P, q, offset, A, b, senses


if __name__ == "__main__":
    # Получение данных для создания модели MILP
    vars_type, lower_bounds, upper_bounds, count, P, q, offset, A, b, senses = get_Data()

    # Инициализация модели
    model = ap.MatrixModel()
    # Добавление переменных
    model.add_variables(
        count=count, lower_bounds=lower_bounds, upper_bounds=upper_bounds, var_types=vars_type
    )
    # Добавление целевой функции
    model.set_objective(P=P, q=q, sense=ap.objective_sense.minimize)
    # Добавление константы к целевой функции
    model.set_objective_offset(offset)
    # Добавление ограничений
    model.add_constraints(A=A, senses=senses, b=b)

    # Сохранение модели в формате LP
    model.write_lp('MIQP.lp')
    # Сохранение модели в формате MPS
    model.write_mps('MIQP.mps')

    # Установление параметров расчета
    model.set_dbl_param('time_limit', 10)
    model.set_bool_param('log_to_console', False)

    # Запуск оптимизации
    result = model.solve()

    # Вывод получившегося результата
    if result.get_solve_result() == ap.solve_result.success:
        if (result.get_solution_status() == ap.solution_status.optimal or \
            result.get_solution_status() == ap.solution_status.feasible):
            print(f'Status problem: {result.get_solve_result()}')
            print(f'Objective value: {result.get_objective_value()}')
            print(f'GAP: {result.get_relative_gap()}')
            print(f'TIME: {result.get_solve_time()} s')
            print(f'Solution vector: {result.get_solution_vector()}')

            # Сохранение найденного решения
            result.write_solution('result.sol')
        else:
            print(f'Status problem: {result.get_solution_status()}')
    else:
        print(f'Status process: {result.get_solve_result()}')

Пример 3

Постановка задачи

Переменные

  • \(x_1 \in \mathbb{B}\) - бинарная переменная,

  • \(x_2 \in \mathbb{Z}\) - целочисленная переменная,

  • \(x_3, x_4 \in \mathbb{R}\) - непрерывные переменные.

Целевая функция

\[ x_1^2 + 2x_2^2 + x_3^2 + 2.5x_1x_4 + x_1 + 3x_3 \ \to \max\limits_{x_1, x_2, x_3, x_4}\]

Ограничения

\[ x_1^2 + x_2^2 + x_3^2 \leqslant 20 \]
\[ x_1 + x_2 + x_3 + x_4 \leqslant 10 \]
\[ x_1 + 2x_2 + 3x_3 \geqslant 5 \]

Решение

import arhiplexpy as ap
import numpy as np
from scipy.sparse import csr_matrix

def get_Data():
    """
    :return: Возвращает данные для создания модели
    """

    vars_type = np.array(
        [ap.variable_type.binary, ap.variable_type.integer, ap.variable_type.continuous, ap.variable_type.continuous]
        , np.int32
    )
    lower_bounds = np.array([0, -ap.kArhiplexInf, -ap.kArhiplexInf, -ap.kArhiplexInf])
    count = 4

    # Диагональ должна быть удвоена относительно исходной матрицы!
    val = np.array([2, 4, 2, 2.5, 2.5])
    row = np.array([0, 1, 2, 0, 3])
    col = np.array([0, 1, 2, 3, 0])
    P = csr_matrix((val, (row, col)), shape=(count, count))
    q = np.array([1, 0, 3, 0])

    A = np.array([
        [1, 1, 1, 1]
        , [1, 2, 3, 0]
    ])
    b = np.array([10, 5])
    senses = np.array(
        [
            ap.constraint_sense.less_equal
            , ap.constraint_sense.greater_equal
        ]
        , np.int32
    )
    val = np.array([2, 2, 2])
    row = np.array([0, 1, 2])
    col = np.array([0, 1, 2])
    P_constr= csr_matrix((val, (row, col)), shape=(count, count))

    return vars_type, lower_bounds, count, P, q, A, b, senses, P_constr


if __name__ == "__main__":
    # Получение данных для создания модели MILP
    vars_type, lower_bounds, count, P, q, A, b, senses, P_constr = get_Data()

    # Инициализация модели
    model = ap.MatrixModel()
    # Добавление переменных
    model.add_variables(
        count=count, lower_bounds=lower_bounds, var_types=vars_type
    )
    # Добавление целевой функции
    model.set_objective(P=P, q=q, sense=ap.objective_sense.maximize)
    # Добавление квадратичного ограничения
    model.add_constraint(P=P_constr, q=None, sense=ap.constraint_sense.less_equal, r=20)
    # Добавление линейных ограничений
    model.add_constraints(A=A, senses=senses, b=b)

    # Сохранение модели в формате LP
    model.write_lp('MIQP.lp')
    # Сохранение модели в формате MPS
    model.write_mps('MIQP.mps')

    # Установление параметров расчета
    model.set_dbl_param('time_limit', 10)
    model.set_bool_param('log_to_console', False)

    # Запуск оптимизации
    result = model.solve()

    # Вывод получившегося результата
    if result.get_solve_result() == ap.solve_result.success:
        if (result.get_solution_status() == ap.solution_status.optimal or \
            result.get_solution_status() == ap.solution_status.feasible):
            print(f'Status problem: {result.get_solve_result()}')
            print(f'Objective value: {result.get_objective_value()}')
            print(f'GAP: {result.get_relative_gap()}')
            print(f'TIME: {result.get_solve_time()} s')
            print(f'Solution vector: {result.get_solution_vector()}')

            # Сохранение найденного решения
            result.write_solution('result.sol')
        else:
            print(f'Status problem: {result.get_solution_status()}')
    else:
        print(f'Status process: {result.get_solve_result()}')