I am currently working on the optimization solution for some industry tasks. I currently use Pyomo and usually write algorithms in Spyder. Due to certain limitations, I have to do it on Windows x64 PCs. I am struggling with Couenne solver. I was able to install the package, my goal was to check how it works. But when I tried to test it on a mixed-integer nonlinear problem with multiply variables and constraints, I got totally different output rather than in the provided example.
For the test I've found one non-linear mixed-integer problem solution tutorial, that is shown here: https://www.youtube.com/watch?v=80M-x1ksiNs&list=PUYeCtS8JUZATtqWsZNKgjVg&index=18
Script (on video) and reference output of the solution: https://www.youtube.com/watch?v=Dl_YLkiRxoc&list=PUYeCtS8JUZATtqWsZNKgjVg&index=22
The problem code was copied from provided instructions: (Maybe I copied it wrong, but I checked it few times)
from pyomo.environ import *
from pyomo.opt import SolverFactory
import numpy as np
# Problem parameters
input_torque = 1400.53333333
fwd_ratio=4.3
rev_ratio=9.1
# Initial diameter guess value
D_guess=np.array([4.5,11.5,6.0,10.0,4.5,2.5,2.25,9.0])
# Initial guess of the teeth number
N_guess=np.array([18,46,24,40,18,20,18,72])
# Initial pitch diameter guess
P_guess =np.array([4,5,18])
# Initial face width
b_guess=np.array([1,1,1])
# Initial central distance guess
C_guess= 5
# Index for number of teeth and the diameter
idx=np.arange(0,8)
# Index for pitch and face width
idx1=np.arange(0,3)
D_dict={}
N_dict={}
P_dict={}
b_dict={}
for i in idx:
D_dict[i]=D_guess[i]
N_dict[i]=N_guess[i]
for i in idx1:
P_dict[i]=P_guess[i]
b_dict[i]=b_guess[i]
# Creating the model object. ConcreteModel instance is appropriate
# as all of the data needed to complete the model is avaiable.
model = ConcreteModel()
# Declare the model variables
# Using Index to access individual gear diameter, bound is the domain
model.D = Var([1,2,3,4,5,6,7,8],initialize=D_dict, within=PositiveReals, bounds=(1.5,30))
model.N = Var([1,2,3,4,5,6,7,8],initialize=N_dict, within=PositiveIntegers, bounds=(18,300))
model.P = Var([1,2,3],initialize=P_dict, within=PositiveIntegers, bounds=(4,80))
model.b = Var([1,2,3],initialize=b_dict, within=PositiveReals)
model.C = Var(initialize=C_guess, within=PositiveReals)
# Declare the objective function
# Represent the total error
def obj_fun(model):
# Defining the residuals
# Error between the model (via diameter) and goal ratio
R1=abs((model.D[2]/model.D[1])*(model.D[4]/model.D[3])-fwd_ratio)
R2=abs((model.D[2]/model.D[5])*(model.D[7]/model.D[6])*(model.D[8]/model.D[7])-rev_ratio)
# Error between the model (via number of teeth) and goal ratio
R3=abs((model.N[2]/model.N[1])*(model.N[4]/model.N[3])-fwd_ratio)
R4=abs((model.N[2]/model.N[5])*(model.N[7]/model.N[6])*(model.N[8]/model.N[7])-rev_ratio)
return R1+R2+R3+R4
# Declare constraints functions
# Diametral pitch constraints
def con_fun1(model):
return model.N[1]/model.D[1] == model.P[1]
def con_fun2(model):
return model.N[2]/model.D[2] == model.P[1]
def con_fun3(model):
return model.N[3]/model.D[3] == model.P[2]
def con_fun4(model):
return model.N[4]/model.D[4] == model.P[2]
def con_fun5(model):
return model.N[5]/model.D[5] == model.P[1]
def con_fun6(model):
return model.N[6]/model.D[6] == model.P[3]
def con_fun7(model):
return model.N[7]/model.D[7] == model.P[3]
def con_fun8(model):
return model.N[8]/model.D[8] == model.P[3]
# Central Distance constraints
def con_fun9(model):
return (model.D[1]+model.D[2])/2 == model.C
def con_fun10(model):
return (model.D[3]+model.D[4])/2 == model.C
def con_fun11(model):
return (model.D[5]+model.D[2])/2 == model.C
def con_fun12(model):
return (model.D[6]+model.D[8])/2 + model.D[7] == model.C
# Declare the pitch-face width constraints
def con_fun13(model):
return 9.0/model.P[1] <= model.b[1]
def con_fun14(model):
return model.b[1] <= 14.0/model.P[1]
def con_fun15(model):
return 9.0/model.P[2] <= model.b[2]
def con_fun16(model):
return model.b[2] <= 14.0/model.P[2]
def con_fun17(model):
return 9.0/model.P[3] <= model.b[3]
def con_fun18(model):
return model.b[3] <= 14.0/model.P[3]
# Declare Gear Tooth Constraint
# Fundamentals of Machine Component Design (Simplified)
def con_fun19(model):
d=model.D[1]
n=model.N[1]
p=model.P[1]
b=model.b[1]
ratio = 1.0
Ft = (input_torque*ratio)/(d/2)
pressure_angle = np.deg2rad(20)
# J - Geometry factor
J = 0.3
# Kv - velocity factor
Kv= 2.5
# Km - mounting factor
Km = 1.4
# Cg - gradient factor
Cg = 0.85
Bn = 388 # AISI 4140 Steel
Kt = 1
Sn = Bn/4
Cl = 1.0
Kr = 0.753
Cs = 0.634
Kms = 1.0
Ko = 1.0
gear_stress = (((Ft*p)/(b*J))*Kv*Ko*Km)/1000
gear_strength = Sn*Cl*Cg*Cs*Kr*Kt*Kms
return gear_stress <= gear_strength
def con_fun20(model):
d=model.D[2]
n=model.N[2]
p=model.P[1]
b=model.b[1]
ratio = model.D[2]/model.D[1]
Ft = (input_torque*ratio)/(d/2)
pressure_angle = np.deg2rad(20)
# J - Geometry factor
J = 0.3
# Kv - velocity factor
Kv= 2.5
# Km - mounting factor
Km = 1.4
# Cg - gradient factor
Cg = 0.85
Bn = 388 # AISI 4140 Steel
Kt = 1
Sn = Bn/4
Cl = 1.0
Kr = 0.753
Cs = 0.634
Kms = 1.0
Ko = 1.0
gear_stress = (((Ft*p)/(b*J))*Kv*Ko*Km)/1000
gear_strength = Sn*Cl*Cg*Cs*Kr*Kt*Kms
return gear_stress <= gear_strength
def con_fun21(model):
d=model.D[3]
n=model.N[3]
p=model.P[2]
b=model.b[2]
ratio = model.D[2]/model.D[1]
Ft = (input_torque*ratio)/(d/2)
pressure_angle = np.deg2rad(20)
# J - Geometry factor
J = 0.3
# Kv - velocity factor
Kv= 2.5
# Km - mounting factor
Km = 1.4
# Cg - gradient factor
Cg = 0.85
Bn = 388 # AISI 4140 Steel
Kt = 1
Sn = Bn/4
Cl = 1.0
Kr = 0.753
Cs = 0.634
Kms = 1.0
Ko = 1.0
gear_stress = (((Ft*p)/(b*J))*Kv*Ko*Km)/1000
gear_strength = Sn*Cl*Cg*Cs*Kr*Kt*Kms
return gear_stress <= gear_strength
def con_fun22(model):
d=model.D[4]
n=model.N[4]
p=model.P[2]
b=model.b[2]
ratio = (model.D[2]/model.D[1])*(model.D[4]/model.D[3])
Ft = (input_torque*ratio)/(d/2)
pressure_angle = np.deg2rad(20)
# J - Geometry factor
J = 0.3
# Kv - velocity factor
Kv= 2.5
# Km - mounting factor
Km = 1.4
# Cg - gradient factor
Cg = 0.85
Bn = 388 # AISI 4140 Steel
Kt = 1
Sn = Bn/4
Cl = 1.0
Kr = 0.753
Cs = 0.634
Kms = 1.0
Ko = 1.0
gear_stress = (((Ft*p)/(b*J))*Kv*Ko*Km)/1000
gear_strength = Sn*Cl*Cg*Cs*Kr*Kt*Kms
return gear_stress <= gear_strength
def con_fun23(model):
d=model.D[5]
n=model.N[5]
p=model.P[1]
b=model.b[1]
ratio = (model.D[2]/model.D[5])
Ft = (input_torque*ratio)/(d/2)
pressure_angle = np.deg2rad(20)
# J - Geometry factor
J = 0.3
# Kv - velocity factor
Kv= 2.5
# Km - mounting factor
Km = 1.4
# Cg - gradient factor
Cg = 0.85
Bn = 388 # AISI 4140 Steel
Kt = 1
Sn = Bn/4
Cl = 1.0
Kr = 0.753
Cs = 0.634
Kms = 1.0
Ko = 1.0
gear_stress = (((Ft*p)/(b*J))*Kv*Ko*Km)/1000
gear_strength = Sn*Cl*Cg*Cs*Kr*Kt*Kms
return gear_stress <= gear_strength
def con_fun24(model):
d=model.D[6]
n=model.N[6]
p=model.P[3]
b=model.b[3]
ratio = (model.D[2]/model.D[5])
Ft = (input_torque*ratio)/(d/2)
pressure_angle = np.deg2rad(20)
# J - Geometry factor
J = 0.3
# Kv - velocity factor
Kv= 2.5
# Km - mounting factor
Km = 1.4
# Cg - gradient factor
Cg = 0.85
Bn = 388 # AISI 4140 Steel
Kt = 1
Sn = Bn/4
Cl = 1.0
Kr = 0.753
Cs = 0.634
Kms = 1.0
Ko = 1.0
gear_stress = (((Ft*p)/(b*J))*Kv*Ko*Km)/1000
gear_strength = Sn*Cl*Cg*Cs*Kr*Kt*Kms
return gear_stress <= gear_strength
def con_fun25(model):
d=model.D[7]
n=model.N[7]
p=model.P[3]
b=model.b[3]
ratio = (model.D[2]/model.D[5])*(model.D[7]/model.D[6])
Ft = (input_torque*ratio)/(d/2)
pressure_angle = np.deg2rad(20)
# J - Geometry factor
J = 0.3
# Kv - velocity factor
Kv= 2.5
# Km - mounting factor
Km = 1.4
# Cg - gradient factor
Cg = 0.85
Bn = 388 # AISI 4140 Steel
Kt = 1
Sn = Bn/4
Cl = 1.0
Kr = 0.753
Cs = 0.634
Kms = 1.0
Ko = 1.0
gear_stress = (((Ft*p)/(b*J))*Kv*Ko*Km)/1000
gear_strength = Sn*Cl*Cg*Cs*Kr*Kt*Kms
return gear_stress <= gear_strength
def con_fun26(model):
d=model.D[8]
n=model.N[8]
p=model.P[3]
b=model.b[3]
ratio = (model.D[2]/model.D[5])*(model.D[8]/model.D[6])
Ft = (input_torque*ratio)/(d/2)
pressure_angle = np.deg2rad(20)
# J - Geometry factor
J = 0.3
# Kv - velocity factor
Kv= 2.5
# Km - mounting factor
Km = 1.4
# Cg - gradient factor
Cg = 0.85
Bn = 388 # AISI 4140 Steel
Kt = 1
Sn = Bn/4
Cl = 1.0
Kr = 0.753
Cs = 0.634
Kms = 1.0
Ko = 1.0
gear_stress = (((Ft*p)/(b*J))*Kv*Ko*Km)/1000
gear_strength = Sn*Cl*Cg*Cs*Kr*Kt*Kms
return gear_stress <= gear_strength
# Assigning model constraints
# Rule - return True or False
model.c1 = Constraint(rule=con_fun1)
model.c2 = Constraint(rule=con_fun2)
model.c3 = Constraint(rule=con_fun3)
model.c4 = Constraint(rule=con_fun4)
model.c5 = Constraint(rule=con_fun5)
model.c6 = Constraint(rule=con_fun6)
model.c7 = Constraint(rule=con_fun7)
model.c8 = Constraint(rule=con_fun8)
model.c9 = Constraint(rule=con_fun9)
model.c10 = Constraint(rule=con_fun10)
model.c11 = Constraint(rule=con_fun11)
model.c12 = Constraint(rule=con_fun12)
model.c13 = Constraint(rule=con_fun13)
model.c14 = Constraint(rule=con_fun14)
model.c15 = Constraint(rule=con_fun15)
model.c16 = Constraint(rule=con_fun16)
model.c17 = Constraint(rule=con_fun17)
model.c18 = Constraint(rule=con_fun18)
model.c19 = Constraint(rule=con_fun19)
model.c20 = Constraint(rule=con_fun20)
model.c21 = Constraint(rule=con_fun21)
model.c22 = Constraint(rule=con_fun22)
model.c23 = Constraint(rule=con_fun23)
model.c24 = Constraint(rule=con_fun24)
model.c25 = Constraint(rule=con_fun25)
model.c26 = Constraint(rule=con_fun26)
# Creating objective function for the model
model.obj = Objective(rule=obj_fun, sense=1)
# Solve the problem
opt = SolverFactory('couenne')
instance = model.create_instance()
results = opt.solve(instance)
instance.solutions.store_to(results)
print(results)
The reference solution was obtained in the Linux system. The solver did it in 30 minutes and was stopped due to the manual initialization of "time out". In my case, the "optimal" solution is always found within 20-30 seconds, but the results are obviously far away from optimum (objective are too high) and integer parameters were outputted as a float.
My output is following:
Problem:
- Lower bound: -inf
Upper bound: inf
Number of objectives: 1
Number of constraints: 0
Number of variables: 23
Sense: unknown
Solver:
- Status: ok
Message: couenne\x3a Optimal
Termination condition: optimal
Id: 3
Error rc: 0
Time: 24.022721529006958
Solution:
- number of solutions: 1
number of solutions displayed: 1
- Gap: None
Status: optimal
Message: couenne\x3a Optimal
Objective:
obj:
Value: 5.3362225559607985
Variable:
C:
Value: 24
D[1]:
Value: 19.75
D[2]:
Value: 28.250000000000004
D[3]:
Value: 22.125000000000007
D[4]:
Value: 25.874999999999993
D[5]:
Value: 19.749999999999996
D[6]:
Value: 4
D[7]:
Value: 9.333333333333332
D[8]:
Value: 25.333333333333336
N[1]:
Value: 79
N[2]:
Value: 113.00000000000001
N[3]:
Value: 177.00000000000006
N[4]:
Value: 206.99999999999997
N[5]:
Value: 78.99999999999997
N[6]:
Value: 24
N[7]:
Value: 56
N[8]:
Value: 152.00000000000006
P[1]:
Value: 4
P[2]:
Value: 8
P[3]:
Value: 6
b[1]:
Value: 3.5
b[2]:
Value: 1.5532917757131675
b[3]:
Value: 1.7813036269011924
Constraint: No values
Can you guess the source of the problem?