Solving linear programs using Pulp but with nonlinear constraint

Question

I am currently sovling an optimization problem using Pulp (I have to using pulp). However, the problem requires the answer of the variables to be perfect squares. Could anyone please help me how to approach with this constraint?

I tried to replace my variable A with B^2 but clearly pulp only supports linear equations. What other ways can I try?

Answer 1

Here is a possible approach. You can use a binary indicator value to indicate which value out of a set of possible values to select. In your case, the set of possible values is the set of perfect squares, but it could be anything.

The complications with this are:

1. You are making a Mixed Integer program, which might be tougher to solve depending on the scope of your problem.
2. You need to do the summation "thing" that I'm showing in the constraint on x anytime you want to capture the value of the squared variable because you need to multiply the binary selection variable by the value to capture it.

aside: you probably do not need to modify the CBC_CMD like I'm doing here... I just need to point to a specific install

Code

# perfect square selector

import pulp

val_index = list(range(10))
square_vals = {t : t**2 for t in val_index}

prob = pulp.LpProblem('picker', sense=pulp.LpMaximize)

x = {t: pulp.LpVariable(f'x_{t}', cat=pulp.LpBinary) for t in val_index}
y = {t: pulp.LpVariable(f'y_{t}', cat=pulp.LpBinary) for t in val_index}

# OBJ:  maximize x^2 + y^2
prob += sum(square_vals[i] * x[i] + square_vals[i] * y[i] for i in val_index)

# CONSTRAINTS

# must select one of the values for x, y
prob += pulp.lpSum(x) == 1
prob += pulp.lpSum(y) == 1

# x^2 <= 50
prob += sum(square_vals[i] * x[i] for i in val_index) <= 50

print(prob)

cbc_path = '/opt/homebrew/opt/cbc/bin/cbc'
solver = pulp.COIN_CMD(path=cbc_path)

solution = prob.solve(solver)

for i in val_index:
    if x[i].varValue >= 0.1:
        print(f'x^2 = {square_vals[i]}')

for i in val_index:
    if y[i].varValue >= 0.1:
        print(f'y^2 = {square_vals[i]}')

Output (truncated):

Result - Optimal solution found

Objective value:                130.00000000
Enumerated nodes:               0
Total iterations:               0
Time (CPU seconds):             0.00
Time (Wallclock seconds):       0.00

Option for printingOptions changed from normal to all
Total time (CPU seconds):       0.00   (Wallclock seconds):       0.00

x^2 = 49
y^2 = 81