使用 MIP、GUROBI 和 MOSEK 解决圆锥优化问题

Question

我编写了以下代码。我需要使用一些求解器来解决这个问题。我尝试使用 mip，特别是 CBC 求解器并得到了一个解决方案，但这不是最佳解决方案，而是一个可行的解决方案。我尝试了各种其他求解器，如 gurobi、mosek、cvxopt、ampl 等，但它们无法解决这个问题，或者由于我对 python 不太熟悉，所以出现了一些复杂情况。请假设所有变量都已定义。

1. 我也想知道为什么当CBC解决这个问题时我没有得到最优值。
2. 我想讨论如何使用求解器来解决此类问题。

如果您需要更多信息，请告诉我。请帮助我。

import numpy as np

# System details
dim = np.shape(A)
rows, col = dim
nv = col  # Total number of variables
dof = nv - rows  # Degree of freedom for the system


def cmatrixcalculate(P_var, A, nv, dof):
    Var = list(range(0, nv))
    Svar = np.array(list(set(Var) - set(P_var)))
    #    print(Var, P_var, Svar)
    Ap = A[:, P_var]
    #    print(Ap)
    As = A[:, Svar]
    #    print(As)
    B = (-1) * np.linalg.inv(As).dot(Ap)
    I = np.eye(dof)

    # Creating the C matrix
    k = 0
    l = 0
    C = np.zeros((nv, B.shape[1]))  # Assuming B has the correct number of columns
    for i in range(nv):
        if i not in P_var:
            k += 1
            C[i, :] = B[k - 1, :]
        else:
            l += 1
            C[i, :] = I[l - 1, :]
    return C


C = cmatrixcalculate(P_var, A, nv, dof)  # Process matrix for the system

inv_sigma = np.zeros(nv)

for ii in range(0, nv):
    inv_sigma[ii] = 1 / sigma_sq[ii]

SIG = np.diag(inv_sigma)

P = np.ones((1, nv))
S = np.linalg.cholesky(C.T @ SIG @ C)
Sigma_g = np.linalg.inv(C.T @ SIG @ C)

from mip import Model, MINIMIZE, CBC, BINARY, OptimizationStatus

# Define the MIP model
model = Model(sense=MINIMIZE,
              solver_name=CBC)  ##CBC is  COIN branch and cut, COIN dedicated for solving optimization problem in industry and academia

# Customize CBC solver parameters
model.emphasis = 1  # 0: default, 1: feasible solutions, 2: proven optimal solution
model.time_limit = 3000000000  # Maximum time limit for solving in seconds

# Define binary variables
x = [model.add_var(var_type=BINARY) for _ in range(nv)]

# Define Q matrix
Q = np.diag(SIG @ x)

CQ = C.T @ Q @ C  # + 1e-6 * np.eye(C.shape[1])

# Define the objective function
model.objective = 0.5 * (dof + dof * np.log(2 * np.pi) - np.linalg.slogdet(CQ.astype(float))[1])

# Add constraint
model += sum(x) == dof
for i in range(dof):
    for j in range(dof):
        model += CQ[i, j] >= 0.1 * np.eye(dof)[i, j]

# Optimize the model
result = model.optimize(max_seconds=1073741824, max_nodes=1073741824, max_solutions=1073741824,
                        max_seconds_same_incumbent=1073741824, max_nodes_same_incumbent=1073741824)

# Create an array to store variable values
variable_values_array = np.zeros(nv)


# Extract values from model.vars and store them in the array
for i, var in enumerate(model.vars):
    variable_values_array[i] = var.x

Q1 = np.diag(SIG @ x)

CQ1 = C.T @ Q @ C

objective_value = 0.5 * (dof + dof * np.log(2 * np.pi) - np.linalg.slogdet(CQ.astype(float))[1])

if result == OptimizationStatus.OPTIMAL:
    print('Optimal function value =  {} '.format(objective_value))
elif result == OptimizationStatus.FEASIBLE:
    print('sol.cost {} found, best possible: {}'.format(model.objective_value, model.objective_bound))
elif result == OptimizationStatus.NO_SOLUTION_FOUND:
    print('no feasible solution found, lower bound is: {}'.format(model.objective_bound))

if result == OptimizationStatus.OPTIMAL or result == OptimizationStatus.FEASIBLE:
    print('Optimal Sensor Network:')
    print(variable_values_array)

Answer 1

即使是小型 MIP 也可能需要花费大量时间来解决。

您的问题很可能属于此类问题。你已经在很大程度上验证了这个结论，因为几个好的优化器都无法解决你的问题。

我建议您发布其中一个优化器的日志输出，因为这通常会提供有关您的问题的宝贵见解。

顺便说一句，您如何得出 CBC 没有报告最佳解决方案的结论？另外，您确定该解决方案不满足 CBC 采用的（可能很弱）停止标准吗？

PS：我记得传感器网络有一些文献。