CODEWITHSUNDEEP
pythonoptimizationmosek
Arjun
Arjun

Reputation: 1

Solving a conic optimization problem using MIP, GUROBI and MOSEK

I have the following code that I wrote. I need to solve this using some solver. I tried using mip, specifically the CBC solver and got a solution, but it wasn't the optimal solution, but a feasible one. I tried various other solvers like gurobi, mosek, cvxopt, ampl etc., but they are not able to solve this or there were complications due to me being new to python in general. The optimla value of the solution to this problem is 13.04.

  1. I would also like to know why I did not get the optimal value when this was solved by CBC.
  2. I would like to discuss how to use solvers for problems like this in general.

Please let me know if you need any more information. Kindly help me.

import numpy as np

A = np.array([[1, 1, -1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
     [0, 0, 0, 0, -1, -1, 1, 1, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
     [-1, 0, 0, 0, 1, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
     [0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
     [0, 0, 1, 0, 0, 0, 0, 0, 0, 0, -1, 0, 1, -1, -1, -1, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
     [0, -1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
     [0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, -1, -1, -1, 0, 0, 0, 0, 0, 0, 0],
     [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, -1, 0, 0, 0, 1, -1, -1, 0, 0, 0, 0],
     [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, -1, 0, 0, -1, 0, 0, 0],
     [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, -1, 1, 0],
     [0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 1]])

sigma_sq = 0.025 * np.array([0.86, 1.00, 111.82, 109.95, 53.27, 112.27, 2.32, 164.05, 0.86, 52.41, 14.86, 67.27, 111.27, 91.86, 60.00, 23.64, 32.73, 16.23, 7.95, 10.50, 87.27, 5.45, 2.59, 46.64, 85.45, 81.32, 70.77, 72.23])

P_var = np.array([1, 2, 5, 6, 7, 8, 12, 14, 15, 16, 17, 19, 20, 21, 23, 24, 27]) - 1

# 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)

The log output is:

Welcome to the CBC MILP Solver
Version: Trunk
Build Date: Oct 28 2021

Starting solution of the Linear programming relaxation problem using Primal Simplex

Clp0024I Matrix will be packed to eliminate 7689 small elements
Coin0506I Presolve 4 (-286) rows, 4 (-24) columns and 10 (-421) elements
Clp1000I sum of infeasibilities 0 - average 0, 2 fixed columns
Coin0506I Presolve 0 (-4) rows, 0 (-4) columns and 0 (-10) elements
Clp0000I Optimal - objective value 0
Clp0000I Optimal - objective value 0
Coin0511I After Postsolve, objective 0, infeasibilities - dual 0 (0), primal 0 (0)
Clp0000I Optimal - objective value 0
Clp0000I Optimal - objective value 0
Clp0000I Optimal - objective value 0
Coin0511I After Postsolve, objective 0, infeasibilities - dual 0 (0), primal 0 (0)
Clp0032I Optimal objective 0 - 0 iterations time 0.012, Presolve 0.01, Idiot 0.01

Starting MIP optimization
Optimal function value = 19.883518735460385
Optimal Sensor Network: [1. 1. 0. 0. 1. 1. 1. 1. 0. 0. 0. 1. 0. 1. 1. 1. 1. 0. 0. 1. 1. 0. 1. 1. 0. 1. 1. 0.]
Cgl0002I 2 variables fixed
Cgl0003I 0 fixed, 0 tightened bounds, 2 strengthened rows, 0 substitutions
Cgl0003I 0 fixed, 0 tightened bounds, 2 strengthened rows, 0 substitutions
Cgl0003I 0 fixed, 0 tightened bounds, 1 strengthened rows, 0 substitutions
Cgl0003I 0 fixed, 0 tightened bounds, 1 strengthened rows, 0 substitutions
Cgl0003I 0 fixed, 1 tightened bounds, 0 strengthened rows, 0 substitutions
Cgl0004I processed model has 3 rows, 3 columns (3 integer (2 of which binary)) and 9 elements
Coin3009W Conflict graph built in 0.001 seconds, density: 9.524%
Cgl0015I Clique Strengthening extended 0 cliques, 0 were dominated
Cbc0045I Nauty: 8 orbits (1 useful covering 2 variables), 1 generators, group size: 2 - sparse size 56 - took 0.001 seconds
Cbc0038I Initial state - 0 integers unsatisfied sum - 0
Cbc0038I Solution found of 0
Cbc0038I Before mini branch and bound, 3 integers at bound fixed and 0 continuous
Cbc0038I Mini branch and bound did not improve solution (0.02 seconds)
Cbc0038I After 0.02 seconds - Feasibility pump exiting with objective of 0 - took 0.00 seconds
Cbc0012I Integer solution of 0 found by feasibility pump after 0 iterations and 0 nodes (0.02 seconds)
Cbc0001I Search completed - best objective 0, took 0 iterations and 0 nodes (0.02 seconds)
Cbc0035I Maximum depth 0, 0 variables fixed on reduced cost
Total time (CPU seconds): 0.03 (Wallclock seconds): 0.03

Process finished with exit code 0

Upvotes: 0

Views: 204

Answers (1)

ErlingMOSEK
ErlingMOSEK

Reputation: 454

Even small MIPs can take a lot of time to solve.

Most likely your problem belongs to that class of problems. You have verified that conclusion to a large extent since several good optimizers cannot solve your problem.

I suggest you post the log output for one of the optimizers because that typically provides valuable insights about your problem.

Btw how do you conclude CBC does not report an optimal solution? Also are you sure the solution does not satisfy the (perhaps weak) stopping criterion CBC employ?

PS: I recall sensor network has some literature.

Upvotes: 0

Related Questions