Gekko discrete optimization indexing problem

Question

I have a question regarding this problem. I want to let model search for optimal value of loader, loader_size, truck, truck_size to solve two constraints of total cost and total duration. However, i faced with a problem regarding the cost. Since truck size and truck, loader size and loader cost cost must be correlated. For example in this code, if the truck size = 20 , truck cost must equal to = 62500 because they are dependent to each other and it is the same way for loader size and loader cost. Please kindly help me with this problem. Thank you.

import numpy as np
from gekko import GEKKO
import random
m = GEKKO()
m.options.SOLVER = 3


Total_Cost = m.Const(15000000)
Total_Duration = m.Const(200)
Q = 10000
k = 1
f = 0.87
E = 0.7
s_load = 7
s_empty = 8
dist = 0.3
dens = 1200
cm = 0.98

#VARIABLE
Loader = m.Var(value = 1, lb=1, ub=2,integer = True)
Truck = m.Var(value = 1, lb=1, ub=5,integer = True)
Loader_Size = m.sos1([2.5, 1.7,1])      
Loader_Cost = m.sos1([21000, 19000,17000])
Truck_Size = m.sos1([20,25.5,15])
Truck_Cost = m.sos1([62500,80000,57000]) 
                             
#DEPENDENT VARIABLE

Travel_T = m.Intermediate((dist/s_load+dist/s_empty)*60)
Unit_Load_Q = m.Intermediate(Truck_Size/dens*1000)
Unit_Load_T = m.Intermediate(Unit_Load_Q/(Loader_Size*k*f*E/cm))
Cycle_Num = m.Intermediate(Q/Unit_Load_Q)
oper_prod = m.min3(Loader/Unit_Load_T,Truck/(Unit_Load_T+Travel_T))
#CONSTRAINTS

a = m.Intermediate(Cycle_Num/oper_prod/60)
b = m.Intermediate((Loader*Loader_Cost)*a)
c = m.Intermediate((Truck*Truck_Cost)*a)
d = m.Intermediate(b+c)


m.Equation(a < Total_Duration)
m.Equation(a > 0)
m.Equation(d < Total_Cost)
m.Equation(d > 0)

#SIMULATION
m.Minimize(d)
m.Minimize(a)
m.solve(disp = False)

print (Loader.value, Loader_Size.value, Loader_Cost.value, Truck.value, Truck_Size.value, Truck_Cost.value)

John Hedengren · Answer 1 · 2022-10-20T13:19:34.287

Use a cubic spline (cspline) to lookup the values of size and cost. They are interpolations between the points but exactly equal the size and cost at the integer values. Here is a script that solves successfully:

import numpy as np
from gekko import GEKKO
import random
m = GEKKO()

Total_Cost = m.Const(15000000)
Total_Duration = m.Const(200)
a = m.Var(lb=0,ub=Total_Duration)
d = m.Var(lb=0,ub=Total_Cost)
Q = 10000
k = 1
f = 0.87
E = 0.7
s_load = 7
s_empty = 8
dist = 0.3
dens = 1200
cm = 0.98

#VARIABLE
Loader = m.Var(value=1, lb=1, ub=5,integer = True)
Loader_Size=m.Var(1.7)
Loader_Cost=m.Var(19000)
Ldr = np.array([1,2,3,4,5])
Ldrs= np.array([2.5,1.7,1,0.5,0.25])
Ldrc= np.array([21000,19000,17000,13000,9000])
m.cspline(Loader,Loader_Size,Ldr,Ldrs,True)
m.cspline(Loader,Loader_Cost,Ldr,Ldrc,True)
#Loader_Size = m.sos1([2.5, 1.7,1])      
#Loader_Cost = m.sos1([21000, 19000,17000])

Truck = m.Var(value=1, lb=1, ub=4,integer = True)
Truck_Size=m.Var(15)
Truck_Cost=m.Var(57000)
Tr = np.array([1,2,3,4])
Trs= np.array([25.5,20,15,10])
Trc= np.array([80000,62500,57000,32000])
m.cspline(Truck,Truck_Size,Tr,Trs,True)
m.cspline(Truck,Truck_Cost,Tr,Trc,True)
#Truck_Size = m.sos1([20,25.5,15,10])
#Truck_Cost = m.sos1([62500,80000,57000,32000]) 
                             
#DEPENDENT VARIABLE

Travel_T = m.Intermediate((dist/s_load+dist/s_empty)*60)
Unit_Load_Q = m.Intermediate(Truck_Size/dens*1000)
Unit_Load_T = m.Intermediate(Unit_Load_Q/(Loader_Size*k*f*E/cm))
Cycle_Num = m.Intermediate(Q/Unit_Load_Q)
oper_prod = m.min3(Loader/Unit_Load_T,Truck/(Unit_Load_T+Travel_T))
#CONSTRAINTS

b = m.Intermediate((Loader*Loader_Cost)*a)
c = m.Intermediate((Truck*Truck_Cost)*a)

# Declare as variables and equations
#m.Equation(a < Total_Duration)
#m.Equation(a > 0)
#m.Equation(d < Total_Cost)
#m.Equation(d > 0)
m.Equation(oper_prod*a==Cycle_Num/60)
m.Equation(d==b+c)

#SIMULATION
m.Minimize(d)
m.Minimize(a)

# initialize with IPOPT
m.options.SOLVER = 3
m.solve(disp = True)

# solve MINLP with APOPT
m.options.SOLVER = 1
m.solve(disp = True)

print(Loader.value, Loader_Size.value, Loader_Cost.value,
      Truck.value, Truck_Size.value, Truck_Cost.value)

The solution is:

  70  1.2917685e+07 6.98e-05 7.47e+07  -0.6 1.30e+00    -  1.00e+00 1.99e-01f  2
  71  1.2917684e+07 2.01e-07 7.07e+05  -1.5 4.44e-01    -  1.00e+00 1.00e+00f  1
  72  1.2917684e+07 3.15e-07 6.33e+06  -2.2 3.66e-02    -  1.00e+00 4.63e-01f  2
  73  1.2917684e+07 4.75e-07 5.95e+04  -3.6 2.71e-02    -  1.00e+00 9.91e-01f  1
  74  1.2917684e+07 1.48e-12 6.89e+01  -5.4 2.94e-03    -  1.00e+00 9.99e-01f  1
  75  1.2917684e+07 1.86e-09 1.38e-03  -7.4 2.10e-05    -  1.00e+00 1.00e+00f  1

Number of Iterations....: 75

                                   (scaled)                 (unscaled)
Objective...............:   1.2917684089800559e+07    1.2917684089800559e+07
Dual infeasibility......:   1.3803296452324758e-03    1.3803296452324758e-03
Constraint violation....:   2.4508488805670487e-12    1.8626451492309570e-09
Complementarity.........:   3.9921377093282655e-08    3.9921377093282655e-08
Overall NLP error.......:   1.1157030658549702e-08    1.3803296452324758e-03


Number of objective function evaluations             = 193
Number of objective gradient evaluations             = 72
Number of equality constraint evaluations            = 193
Number of inequality constraint evaluations          = 0
Number of equality constraint Jacobian evaluations   = 77
Number of inequality constraint Jacobian evaluations = 0
Number of Lagrangian Hessian evaluations             = 75
Total CPU secs in IPOPT (w/o function evaluations)   =      0.038
Total CPU secs in NLP function evaluations           =      0.012

EXIT: Optimal Solution Found.
 
 The solution was found.
 
 The final value of the objective function is    12917684.0898006     
 
 ---------------------------------------------------
 Solver         :  IPOPT (v3.12)
 Solution time  :   7.190000000264263E-002 sec
 Objective      :    12917684.0898006     
 Successful solution
 ---------------------------------------------------
 
apm 136.36.4.38_gk_model0 <br><pre> ----------------------------------------------------------------
 APMonitor, Version 1.0.1
 APMonitor Optimization Suite
 ----------------------------------------------------------------
 
 
 --------- APM Model Size ------------
 Each time step contains
   Objects      :            4
   Constants    :            2
   Variables    :           12
   Intermediates:            6
   Connections  :            8
   Equations    :           13
   Residuals    :            7
 
 Number of state variables:             12
 Number of total equations: -            9
 Number of slack variables: -            2
 ---------------------------------------
 Degrees of freedom       :              1
 
 ----------------------------------------------
 Steady State Optimization with APOPT Solver
 ----------------------------------------------
Iter:     1 I:  0 Tm:      0.00 NLPi:   11 Dpth:    0 Lvs:    3 Obj:  1.29E+07 Gap:       NaN
--Integer Solution:   1.47E+07 Lowest Leaf:   1.29E+07 Gap:   1.26E-01
Iter:     2 I:  0 Tm:      0.00 NLPi:    4 Dpth:    1 Lvs:    2 Obj:  1.47E+07 Gap:  1.26E-01
Iter:     3 I:  0 Tm:      0.00 NLPi:   12 Dpth:    1 Lvs:    4 Obj:  1.34E+07 Gap:  1.26E-01
--Integer Solution:   1.47E+07 Lowest Leaf:   1.34E+07 Gap:   9.18E-02
Iter:     4 I:  0 Tm:      0.00 NLPi:    6 Dpth:    1 Lvs:    3 Obj:  1.47E+07 Gap:  9.18E-02
Iter:     5 I:  0 Tm:      0.00 NLPi:    5 Dpth:    2 Lvs:    4 Obj:  1.29E+07 Gap:  9.18E-02
Iter:     6 I: -1 Tm:      0.00 NLPi:    3 Dpth:    3 Lvs:    3 Obj:  1.29E+07 Gap:  9.18E-02
Iter:     7 I: -1 Tm:      0.00 NLPi:    2 Dpth:    3 Lvs:    2 Obj:  1.29E+07 Gap:  9.18E-02
Iter:     8 I: -1 Tm:      0.00 NLPi:    3 Dpth:    2 Lvs:    1 Obj:  1.34E+07 Gap:  9.18E-02
Iter:     9 I: -1 Tm:      0.00 NLPi:    8 Dpth:    2 Lvs:    0 Obj:  1.34E+07 Gap:  9.18E-02
 No additional trial points, returning the best integer solution
 Successful solution
 
 ---------------------------------------------------
 Solver         :  APOPT (v1.0)
 Solution time  :   2.300000000104774E-002 sec
 Objective      :    14654721.8680576     
 Successful solution
 ---------------------------------------------------
 
[1.0] [2.5] [21000.0] [1.0] [25.5] [80000.0]

IPOPT finds an initial non-integer solution and then APOPT finds the integer solution.

It seems like in cspline, we need to add more data point to the array which initially not provided. The code did run but seems like it cant give the most optimal answer at all. I wonder if there is any other way beside using cspline? — kennnn, Oct 20 '22 at 07:01
also it seems like the answer keeps giving the initial value of Loader_Size=m.Var(2.5), Loader_Cost=m.Var(21000) , Truck_Size=m.Var(20) Truck_Cost=m.Var(62500). — kennnn, Oct 20 '22 at 07:18
There are a minimum of 4 data points to use cspline. If you need more to reach 4, just replicate some of the options. I adjusted the initial conditions to show that the values change with the optimization. — John Hedengren, Oct 20 '22 at 13:20

Gekko discrete optimization indexing problem

1 Answers1