I have an optimization problem, that I have to find the the lowest cost of the given motors.
And there is a constraint that, the motor either run or doesn't run. But if it runs, it has to reach the lower limit of its power range
I am going to include my code, to show what I have tried.
from gekko import GEKKO
power_ranges = {
'Motor1': (0.6, 1.1),
'Motor2': (2.1, 6),
'Motor3': (1, 1.94),
'Motor4': (1, 1.94),
}
prices = {
'Motor1': lambda x: (x ** 2)/x*5000,
'Motor2': lambda y: (y ** 1 / 0.45) * 5500,
'Motor3': lambda z: (z * 0.45) * 5100,
'Motor4': lambda a: (a / 0.45) * 5200,
}
model = GEKKO()
# Define decision variables
x = {}
y = {}
for motor in power_ranges:
x[motor] = model.Var(lb=power_ranges[motor][0], ub=power_ranges[motor][1])
y[motor] = model.Var(lb=0, ub=1, integer=True)
# Define objective function
model.Minimize(sum(prices[motor](x[motor]) for motor in power_ranges))
# Define lower and upper bounds constraints
for motor in power_ranges:
lower_bound, upper_bound = power_ranges[motor]
model.Equation(x[motor] <= upper_bound)
# Define the "or" constraint
for motor in power_ranges:
model.Equation(x[motor] >= power_ranges[motor][0] * y[motor])
model.Equation(x[motor] <= power_ranges[motor][1] * y[motor])
# Define power constraint
model.Equation(sum(x[motor] for motor in power_ranges) == 4.7)
# Solve the optimization problem
model.options.SOLVER = 1
model.solve()
# Print the solution
if model.options.APPSTATUS == 1:
print("Optimal solution found:")
for motor in power_ranges:
print(f"{motor}: {round(float(x[motor].value[0]), 2)}")
print(f"Total cost: {round(float(model.options.OBJFCNVAL), 2)}")
else:
print("No optimal solution found.")
And this always gives me an answer where every motor run. Like if run this, i get a solution, but if I change the code to this:
# Define power constraint
model.Equation(sum(x[motor] for motor in power_ranges) == 1)
Or change the power constraint to something which is less tha 4.7 and greater or equal 1 I do not get a solution, but as we can see, if a change the power constraint there are 3 available solution, without price
Thank you in advance.
The equations need to be adjusted to "turn off" the constraint when the motor is not selected when y[motor]=0. The binary switch variable y is the correct way to implement the equation but it is infeasible if not multiplied on the left side of the equation.
# Define the "or" constraint
for motor in power_ranges:
model.Equation(x[motor]* y[motor] >= power_ranges[motor][0] * y[motor])
model.Equation(x[motor]* y[motor] <= power_ranges[motor][1] * y[motor])
# Define power constraint
model.Equation(sum(x[motor]*y[motor] for motor in power_ranges) == 10.98)
One other issue is that the objective function goes to Inf when Motor1 is at zero. A small modification makes the objective go to zero instead of infinity. You can adjust the small denominator term 1e-3 to ensure convergence without changing the original cost function too much. Another solution with m.if3() is shown in the answer to your follow-up question: How to use If...else in lambda with gekko
prices = {
'Motor1': lambda x: (x ** 2)/(x+1e-3)*5000,
'Motor2': lambda y: (y ** 1 / 0.45) * 5500,
'Motor3': lambda z: (z * 0.45) * 5100,
'Motor4': lambda a: (a / 0.45) * 5200,
}
An intermediate variable z[motor] is added to improve the readability of the model. It now gives the optimal solution for power constraint 0, 1...10.98. Above 10.98 or >0 to <1 there is no feasible solution because no combination of motors can provide that power. Here is the complete solution with a modification to the print statement at the bottom:
from gekko import GEKKO
power_ranges = {
'Motor1': (0.6, 1.1),
'Motor2': (2.1, 6),
'Motor3': (1, 1.94),
'Motor4': (1, 1.94),
}
prices = {
'Motor1': lambda x: (x ** 2)/(x+1e-3)*5000,
'Motor2': lambda y: (y ** 1 / 0.45) * 5500,
'Motor3': lambda z: (z * 0.45) * 5100,
'Motor4': lambda a: (a / 0.45) * 5200,
}
model = GEKKO()
# Define decision variables
x = {}; y = {}; z = {}
for motor in power_ranges:
x[motor] = model.Var(lb=power_ranges[motor][0],
ub=power_ranges[motor][1])
y[motor] = model.Var(lb=0, ub=1, integer=True)
z[motor] = model.Intermediate(x[motor]*y[motor])
# Define objective function
model.Minimize(sum(prices[motor](z[motor]) \
for motor in power_ranges))
# Define the "or" constraint
for motor in power_ranges:
model.Equation(z[motor] >= power_ranges[motor][0] * y[motor])
model.Equation(z[motor] <= power_ranges[motor][1] * y[motor])
# Define power constraint
model.Equation(sum(z[motor] for motor in power_ranges) == 4.7)
# Solve the optimization problem
model.options.SOLVER = 1
model.solve()
# Print the solution
if model.options.APPSTATUS == 1:
print("Optimal solution found:")
for motor in power_ranges:
print(f"{motor}: {round(float(z[motor].value[0]), 2)}")
print(f"Total cost: {round(float(model.options.OBJFCNVAL), 2)}")
else:
print("No optimal solution found.")
The optimal solutions with 4.7 for the power demand:
Optimal solution found:
Motor1: 1.1
Motor2: 0.0
Motor3: 1.94
Motor4: 1.66
Total cost: 29129.53