Linear System Constraint Optimization with Integer Solutions

Question

I am interested in solving the following optimization problem and obtaining the lowest positive integer coefficients:

min(x1 + x2 + x3)
s.t. 
2x1 – 2x3 = 0
2x2 – x3 = 0

Since the number of equations do not match the number of variables, I also add the constraint that x1 = 1. Using the documentation from ortools, I adapted the code as follows:

mat = np.array([[2,0,-2],[0,2,-1]])
def create_data_model(mat):
    num_rows, num_vars =mat.shape
    data = {}
    data['constraint_coeffs'] = [
          [2, 0, -2],
          [0, 2, -1],
          [1, 0, 0],
      ]
    data['bounds'] = [0, 0, 1]
    data['obj_coeffs'] = [1, 1, 1]
    data['num_vars'] = num_vars
    offset = np.max((num_rows, num_vars))-np.min((num_rows, num_vars))
    data['num_constraints'] = num_rows+offset
    return data

data=create_data_model(mat)
infinity = solver.infinity()
x = {}
for j in range(data['num_vars']):
    x[j] = solver.IntVar(0, infinity, 'x[%i]' % j)
print('Number of variables =', solver.NumVariables())

for i in range(data['num_constraints']):
    constraint = solver.RowConstraint(0, data['bounds'][i], '')
    for j in range(data['num_vars']):
        constraint.SetCoefficient(x[j], data['constraint_coeffs'][i][j])

print('Number of constraints =', solver.NumConstraints())

objective = solver.Objective()
for j in range(data['num_vars']):
    objective.SetCoefficient(x[j], data['obj_coeffs'][j])
objective.SetMinimization()

if status == pywraplp.Solver.OPTIMAL:
    for j in range(data['num_vars']):
        print(x[j].name(), ' = ', x[j].solution_value())
    print()
    print('Problem solved in %f milliseconds' % solver.wall_time())
    print('Problem solved in %d iterations' % solver.iterations())
    print('Problem solved in %d branch-and-bound nodes' % solver.nodes())
else:
    print('The problem does not have an optimal solution.')

denoms=[x[0].solution_value().as_integer_ratio()[1], x[1].solution_value().as_integer_ratio()[1],
       x[2].solution_value().as_integer_ratio()[1]]
scaling_factor=lcmm(*denoms) #a method to get the least common multiple
coeff=[scaling_factor*x[0].solution_value(), scaling_factor*x[1].solution_value(),
       scaling_factor*x[2].solution_value()]

This should return [1,0.5, 1] as the optimal solution, which becomes [2,1,2] after scaling.

(1) In my implementation, the solver does not find an optimal solution. Why is that?

(2) Is there a simpler way to get to the least positive integer solution?

Answer 1

I see a few issues in your code.

With the line:

x[j] = solver.IntVar(0.0, infinity, 'x[%i]' % j)

You are saying that the optimizer must produce integers. However, the answer to your example was going to give x2 as 0.5, so it will fail to give the optimal solution. Try:

x[j] = solver.NumVar(0.0, infinity, 'x[%i]' % j)

Which allows floats. Then you can multiply it up to ints later.

Another issue, you set each constraint with:

constraint = solver.Constraint(0, data['bounds'][i], '')

Which means the lower bound for this expression is 0, and the upper bound is the bound value in your data structure. However, for your third constraint, this would be like saying 0 <= x1 <= 1 (bounds[2] is 1). However, you wanted to set this as an equality, not a range from 0 to some bound. You can replace the line with

constraint = solver.Constraint(data['bounds'][i], data['bounds'][i], '')

To guarantee that the value is exactly what you want. You can also rename "bounds" to "target" to make it more clear that this isn't an upper or lower bound, but an exact target you want.

Finally, I don't see where you run status = solver.Solve() in your code. It should be right before you check status in that if/else.

Hope this helps!