linear-programminggekkomixed-integer-programming

Gekko Code for Resource Optimization problem


I am working on resource optimization problem and writing a GEKKO code to solve. Problem statement is as follows: Let assume there are 2 workers and 3 tasks. Each worker gets some rewards for the tasks it picks. Each task can be assigned to only one worker and each worker can only pick one tasks. The objective is to maximize the total rewards obtained by both workers.

m = GEKKO()

# decision variables
x11 = m.Var(value=0, lb=0,ub=1,integer=True) # variable if worker 1 picks task 1
x12 = m.Var(value=0, lb=0,ub=1,integer=True) # variable if worker 1 picks task 2
x13 = m.Var(value=0, lb=0,ub=1,integer=True) # variable if worker 1 picks task 3
x21 = m.Var(value=0, lb=0,ub=1,integer=True) # variable if worker 2 picks task 1
x22 = m.Var(value=0, lb=0,ub=1,integer=True) # variable if worker 2 picks task 2
x23 = m.Var(value=0, lb=0,ub=1,integer=True) # variable if worker 2 picks task 3


# rewards constants
S11 = m.Const(value= 2) # reward if worker 1 picks task 1
S12 = m.Const(value= 5) # reward if worker 1 picks task 2
S13 = m.Const(value= 8) # reward if worker 1 picks task 3
S21 = m.Const(value= 1) # reward if worker 2 picks task 1
S22 = m.Const(value= 5) # reward if worker 2 picks task 2
S23 = m.Const(value= 7) # reward if worker 2 picks task 3

# Equations
m.Equation(x11 + x12 + x13 == 1) # constraints for worker 1
m.Equation(x21 + x22 + x23 == 1) # constraints for worker 2

m.Obj(-(x11*S11 + x12*S12 + x13*S13 + x21*S21 + x22*S22 + x23*S23)) # Objective

m.options.IMODE = 3 # Steady state optimization
m.options.SOLVER = 1        # change solver (1=APOPT,3=IPOPT)
m.solve(disp=True)          # solve locally (remote=False)

print('Results')
print('x11: ' + str(x11.value))
print('x12: ' + str(x12.value))
print('x13: ' + str(x13.value))
print('x21: ' + str(x21.value))
print('x22: ' + str(x22.value))
print('x23: ' + str(x23.value))
print('Objective: ' + str(m.options.objfcnval))

upon running this I am getting output as:

    apm 122.172.80.123_gk_model3 <br><pre> ----------------------------------------------------------------
 APMonitor, Version 1.0.1
 APMonitor Optimization Suite
 ----------------------------------------------------------------




 --------- APM Model Size ------------
 Each time step contains
   Objects      :            0
   Constants    :            6
   Variables    :            6
   Intermediates:            0
   Connections  :            0
   Equations    :            3
   Residuals    :            3
 
 Number of state variables:              6
 Number of total equations: -            2

 Number of slack variables: -            0
 ---------------------------------------
 Degrees of freedom       :              4
 
 ----------------------------------------------
 Steady State Optimization with APOPT Solver
 ----------------------------------------------
Iter:     1 I:  0 Tm:      0.00 NLPi:    3 Dpth:    0 Lvs:    0 Obj: -1.50E+01 Gap:  0.00E+00



Successful solution
 
 ---------------------------------------------------
 Solver         :  APOPT (v1.0)
 Solution time  :   1.319999998668209E-002 sec
 Objective      :   -15.0000000000000     
 Successful solution

 ---------------------------------------------------
 
Results
x11: [0.0]
x12: [0.0]
x13: [1.0]
x21: [0.0]
x22: [0.0]
x23: [1.0]
Objective: -15.0

The solution looks correct since I have not added the constraint of each tasks to be assigned to be only to 1 worker. But I add following tasks constraints as :

m.Equation(x11 + x21 == 1) # constraint for task1
m.Equation(x12 + x22 == 1) # constraint for task2
m.Equation(x13 + x23 == 1) # constraint for task3

I am getting following error:

    apm 122.172.80.123_gk_model4 <br><pre> ----------------------------------------------------------------
 APMonitor, Version 1.0.1
 APMonitor Optimization Suite
 ----------------------------------------------------------------
 
 
 --------- APM Model Size ------------
 Each time step contains
   Objects      :            0
   Constants    :            6

   Variables    :            6

   Intermediates:            0
   Connections  :            0
   Equations    :            6
   Residuals    :            6
 
 Number of state variables:              6
 Number of total equations: -            5

 Number of slack variables: -            0
 ---------------------------------------
 Degrees of freedom       :              1
 
 ----------------------------------------------
 Steady State Optimization with APOPT Solver
 ----------------------------------------------
Iter:     1 I: -1 Tm:      0.00 NLPi:    2 Dpth:    0 Lvs:    0 Obj:  0.00E+00 Gap:  



    NaN
 Warning: no more possible trial points and no integer solution
 Maximum iterations
 
 ---------------------------------------------------
 Solver         :  APOPT (v1.0)
 Solution time  :   1.320000001578592E-002 sec
 Objective      :   0.000000000000000E+000
 Unsuccessful with error code            0
 ---------------------------------------------------
 
 Creating file: infeasibilities.txt
 Use command apm_get(server,app,'infeasibilities.txt') to retrieve file
 @error: Solution Not Found
---------------------------------------------------------------------------
Exception                                 Traceback (most recent call last)
Cell In[40], line 14
     12 m.options.IMODE = 3 # Steady state optimization
     13 m.options.SOLVER = 1        # change solver (1=APOPT,3=IPOPT)
---> 14 m.solve(disp=True)          # solve locally (remote=False)
     16 print('Results')
     17 print('x11: ' + str(x11.value))

File ~/anaconda3/lib/python3.11/site-packages/gekko/gekko.py:2185, in GEKKO.solve(self, disp, debug, GUI, **kwargs)
   2183 #print APM error message and die
   2184 if (debug >= 1) and ('@error' in response):
-> 2185     raise Exception(response)
   2187 #load results
   2188 def byte2str(byte):

Exception:  @error: Solution Not Found

Since this is a very trivial example and the solution does exist. Could you please help in understanding what could be going wrong in the code? and also how do we debug such errors?

Thanks


Solution

  • The second set of constraints should not be equality constraints, but inequality constraints. Reproducing on PuLP,

    import pulp
    
    worker_idx = range(1, 3)
    task_idx = range(1, 4)
    
    assigns = pulp.LpVariable.matrix(
        name='assign',
        cat=pulp.LpBinary,
        indices=(worker_idx, task_idx),
    )
    
    rewards = (
        (2, 5, 8),
        (1, 5, 7),
    )
    
    m = pulp.LpProblem(name='resources', sense=pulp.LpMaximize)
    m.setObjective(
        pulp.lpSum(
            pulp.lpDot(worker_assigns, worker_rewards)
            for worker_assigns, worker_rewards in zip(assigns, rewards)
        )
    )
    
    # Each worker can only pick one task
    for worker, worker_assign in zip(worker_idx, assigns):
        m.addConstraint(
            name=f'excl_w{worker}', constraint=pulp.lpSum(worker_assign) == 1,
        )
    
    # Each task can be assigned up to once
    for i, task in enumerate(task_idx):
        m.addConstraint(
            name=f'excl_t{task}', constraint=pulp.lpSum(
                worker_assigns[i]
                for worker_assigns in assigns
            ) <= 1,
        )
    
    print(m)
    m.solve()
    assert m.status == pulp.LpStatusOptimal
    for worker, worker_assigns in zip(worker_idx, assigns):
        task, = (
            task
            for task, assign in zip(task_idx, worker_assigns)
            if assign.value() > 0.5
        )
        print(f'Worker {worker}: task {task}')
    
    resources:
    MAXIMIZE
    2*assign_1_1 + 5*assign_1_2 + 8*assign_1_3 + 1*assign_2_1 + 5*assign_2_2 + 7*assign_2_3 + 0
    SUBJECT TO
    excl_w1: assign_1_1 + assign_1_2 + assign_1_3 = 1
    
    excl_w2: assign_2_1 + assign_2_2 + assign_2_3 = 1
    
    excl_t1: assign_1_1 + assign_2_1 <= 1
    
    excl_t2: assign_1_2 + assign_2_2 <= 1
    
    excl_t3: assign_1_3 + assign_2_3 <= 1
    
    VARIABLES
    0 <= assign_1_1 <= 1 Integer
    0 <= assign_1_2 <= 1 Integer
    0 <= assign_1_3 <= 1 Integer
    0 <= assign_2_1 <= 1 Integer
    0 <= assign_2_2 <= 1 Integer
    0 <= assign_2_3 <= 1 Integer
    
    Welcome to the CBC MILP Solver 
    Version: 2.10.3 
    Build Date: Dec 15 2019 
    
    At line 2 NAME          MODEL
    At line 3 ROWS
    At line 10 COLUMNS
    At line 41 RHS
    At line 47 BOUNDS
    At line 54 ENDATA
    Problem MODEL has 5 rows, 6 columns and 12 elements
    Coin0008I MODEL read with 0 errors
    Option for timeMode changed from cpu to elapsed
    Continuous objective value is 13 - 0.00 seconds
    Cgl0004I processed model has 5 rows, 6 columns (6 integer (6 of which binary)) and 12 elements
    
    
    Result - Optimal solution found
    
    Objective value:                13.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
    
    Worker 1: task 3
    Worker 2: task 2