MILP

=MILP(c, [integrality], [bounds_lower], [bounds_upper], [a_ub], [b_ub], [a_eq], [b_eq])

=MILP({0,-1}, {1,1}, {0,0}, {5,5}, {-1,1;3,2;2,3}, {1;12;12})

=MILP({1,2}, {1,1}, {0,0}, {4,3}, {-1,-1}, {-3})

=MILP({3,2}, {0,1}, , , , , {1,1}, {5}, {0,0}, {5,5})

=MILP({1,1,1}, , , , {1,2,0;0,1,1}, {4;3})

import math
from typing import List, Optional, Union
 
import numpy as np
from scipy.optimize import Bounds, LinearConstraint, milp as scipy_milp
 
 
def milp(
    c: Union[List[List[float]], float],  # Can be 2D list or scalar
    integrality: Optional[Union[List[List[int]], int]] = None,  # Can be 2D list or scalar
    bounds_lower: Optional[Union[List[List[float]], float]] = None,  # Can be 2D list or scalar
    bounds_upper: Optional[Union[List[List[float]], float]] = None,  # Can be 2D list or scalar
    a_ub: Optional[Union[List[List[float]], float]] = None,  # Can be 2D list or scalar
    b_ub: Optional[Union[List[List[float]], float]] = None,  # Can be 2D list or scalar
    a_eq: Optional[Union[List[List[float]], float]] = None,  # Can be 2D list or scalar
    b_eq: Optional[Union[List[List[float]], float]] = None,  # Can be 2D list or scalar
):
    """
    Solve a mixed-integer linear program using scipy.optimize.milp.
    Wrapper for scipy.optimize.milp function (see: https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.milp.html).
 
    Args:
        c: 2D list representing the objective coefficients (row or column vector).
        integrality: 2D list matching the length of c that specifies variable integrality (0=continuous, 1=integer, 2=semi-continuous, 3=semi-integer).
        bounds_lower: 2D list providing lower bounds for each variable.
        bounds_upper: 2D list providing upper bounds for each variable.
        a_ub: 2D list for inequality constraint coefficients (A_ub x <= b_ub).
        b_ub: 2D list for inequality constraint limits.
        a_eq: 2D list for equality constraint coefficients (A_eq x == b_eq).
        b_eq: 2D list for equality constraint limits.
 
    Returns:
        list[list[float]] or str: [[x1, x2, ..., optimal_value]] if successful, otherwise an error message string.
 
    This example function is provided as-is without any representation of accuracy.
    """
    # Validate and flatten objective coefficients
    # Normalize scalar inputs (Excel passes [[x]] as x) to 2D lists
    def normalize_to_2d_list(value):
        if not isinstance(value, list):
            return [[value]]
        return value
 
    c = normalize_to_2d_list(c)
    if integrality is not None:
        integrality = normalize_to_2d_list(integrality)
    if bounds_lower is not None:
        bounds_lower = normalize_to_2d_list(bounds_lower)
    if bounds_upper is not None:
        bounds_upper = normalize_to_2d_list(bounds_upper)
    if a_ub is not None:
        a_ub = normalize_to_2d_list(a_ub)
    if b_ub is not None:
        b_ub = normalize_to_2d_list(b_ub)
    if a_eq is not None:
        a_eq = normalize_to_2d_list(a_eq)
    if b_eq is not None:
        b_eq = normalize_to_2d_list(b_eq)
 
    try:
        c_vec = np.array(c, dtype=float).flatten()
    except Exception:
        return "Invalid input: c must be a 2D list of numeric values."
 
    if c_vec.size == 0:
        return "Invalid input: c must contain at least one coefficient."
 
    n_vars = c_vec.size
 
    # Helper to extract 1D arrays from potential 2D lists
    def _to_array(values, name):
        if values is None:
            return None
        try:
            arr = np.array(values, dtype=float)
        except Exception:
            raise ValueError(f"Invalid input: {name} must contain numeric values.")
        return arr
 
    # Integrality handling
    integrality_arr = None
    if integrality is not None:
        try:
            integrality_arr = np.array(integrality, dtype=int).flatten()
        except Exception:
            return "Invalid input: integrality must be integers."
        if integrality_arr.size != n_vars:
            return "Invalid input: integrality length must match number of variables."
        if not np.all(np.isin(integrality_arr, [0, 1, 2, 3])):
            return "Invalid input: integrality values must be 0 (continuous), 1 (integer), 2 (semi-continuous), or 3 (semi-integer)."
 
    # Bounds handling
    lower_arr = None
    upper_arr = None
    try:
        if bounds_lower is not None:
            lower_arr = _to_array(bounds_lower, "bounds_lower").flatten()
        if bounds_upper is not None:
            upper_arr = _to_array(bounds_upper, "bounds_upper").flatten()
    except ValueError as exc:
        return str(exc)
 
    if lower_arr is not None and lower_arr.size != n_vars:
        return "Invalid input: bounds_lower must provide one value per variable."
    if upper_arr is not None and upper_arr.size != n_vars:
        return "Invalid input: bounds_upper must provide one value per variable."
 
    if lower_arr is not None and upper_arr is not None:
        if np.any(lower_arr > upper_arr):
            return "Invalid input: each lower bound must be less than or equal to the corresponding upper bound."
 
    bounds_obj = None
    if lower_arr is not None or upper_arr is not None:
        # SciPy requires both vectors; if one missing, fill with defaults
        if lower_arr is None:
            lower_arr = np.zeros(n_vars)
        if upper_arr is None:
            upper_arr = np.full(n_vars, math.inf)
        bounds_obj = Bounds(lower_arr, upper_arr)
 
    # Constraint processing helper
    def _build_linear_constraint(matrix, vector, sense):
        if matrix is None and vector is None:
            return None
        if matrix is None or vector is None:
            return f"Invalid input: {sense} constraints require both matrix and vector."
        mat = _to_array(matrix, f"{sense} matrix")
        vec = _to_array(vector, f"{sense} vector").flatten()
        if mat.ndim == 1:
            mat = mat.reshape(1, -1)
        if mat.shape[1] != n_vars:
            return f"Invalid input: {sense} matrix must have {n_vars} columns."
        if vec.size != mat.shape[0]:
            return f"Invalid input: {sense} vector length must equal number of rows in the matrix."
        if sense == "inequality":
            return LinearConstraint(mat, -np.inf, vec)
        return LinearConstraint(mat, vec, vec)
 
    constraints = []
    for sense, matrix, vector in (
        ("inequality", a_ub, b_ub),
        ("equality", a_eq, b_eq),
    ):
        constraint = _build_linear_constraint(matrix, vector, sense)
        if isinstance(constraint, str):
            return constraint
        if constraint is not None:
            constraints.append(constraint)
 
    if len(constraints) == 0:
        constraints = None
 
    try:
        result = scipy_milp(
            c_vec,
            integrality=integrality_arr,
            bounds=bounds_obj,
            constraints=constraints,
        )
    except Exception as exc:
        return f"Error during MILP solving: {exc}"
 
    if not result.success:
        message = result.message if hasattr(result, "message") else "MILP solver did not succeed."
        return f"MILP solving failed: {message}"
 
    solution = result.x
    if solution is None:
        return "MILP solving failed: no solution returned."
 
    try:
        solution_list = list(map(float, solution))
    except Exception:
        return "Error converting MILP solution to floats."
 
    objective = float(result.fun) if result.fun is not None else float(np.dot(c_vec, solution))
 
    return [solution_list + [objective]]